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I We argue that an astrophysical solution to the ultra high energy cosmic ray (UHECR) prob- 

. lem is viable. The detailed study of UHECR energy spectra is performed. The spectral features 

' of extragalactic protons interacting with Cosmic Microwave Background (CMB) are calculated in 

^S) ' model-independent way. Using the power-law generation spectrum oc E~''^ as the only assumption, 

we analyze four features of the proton spectrum: the GZK cutoff, dip, bump and the second dip. We 
^ ' found the dip, induced by electron-positron production on CMB, as the most robust feature, existing 

' in energy range 1 x 10^* — 4 x 10^^ eV. Its shape is stable relative to various phenomena included 

in calculations: discreteness of the source distribution, different modes of UHE proton propagation 
(from rectilinear to diffusive), local overdensity or deficit of the sources, large-scale inhomogeneities 
in the universe and interaction fluctuations. The dip is well confirmed by observations of AGASA, 
I HiRes, Fly's Eye and Yakutsk detectors. With two free parameters (7g and flux normalization 

^ ■ constant) the dip describes about 20 energy bins with x^/d.o.f. « 1 for each experiment. The best 

fit is reached at 7g = 2.7, with the allowed range 2.55 - 2.75. The dip is used for energy calibration 
of the detectors. For each detector independently the energy is shifted by factor A to reach the 
minimum x^- We found AAg = 0.9, Xm = 1.2 and Ava = 0.75 for AGASA, HiRes and Yakutsk 
detectors, respectively. Remarkably, that after this energy shift the fluxes and spectra of all three 
detectors agree perfectly, with discrepancy between AGASA and HiRes at _E > 1 x 10^" eV being 
not statistically significant. The excellent agreement of the dip with observations should be consid- 
ered as confirmation of UHE proton interaction with CMB. The dip has two flattenings. The high 
energy flattening at -E « 1 x 10^^ eV automatically explains ankle, the feature observed in all exper- 
Oh, iments starting from 1980s. The low-energy flattening at iS ~ 1 x 10^* eV provides the transition 

to galactic cosmic rays. This transition is studied quantitatively in this work. Inclusion of primary 
nuclei with fraction more than 20% upsets the agreement of the dip with observations, which we 
^ \ interpret as an indication to acceleration mechanism. We study in detail the formal problems of 

spectra calculations: energy losses (the new detailed calculations are presented) , analytic method of 
spectrum calculations and study of fluctuations with help of kinetic equation. The UHECR sources, 
k>( \ AGN and GRBs, are studied in a model-dependent way, and acceleration is discussed. Based on 

, the agreement of the dip with existing data, we make the robust prediction for the spectrum at 

' 1 X 10^* — 1 X 10^° eV to be measured in the nearest future by Auger detector. We also predict the 

spectral signature of nearby sources, if they are observed by Auger. This paper is long and contains 
many technical details. For those who are interested only in physical content we recommend to read 
Introduction and Conclusions, which are written as autonomous parts of the paper. 

PACS numbers: 98.70.Sa, 96.50.sb, 96.50.sd, 98.54.Cm 
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I. INTRODUCTION 



The systematic study of Ultra High Energy Cosmic Rays (UHECR) started in late 1950s after construction of 
Volcano Ranch (USA) and Moscow University (USSR) arrays. During the next 50 years of research the origin of UHE 
particles, which hit the detectors, was not well understood. At present due to the data of the last generation arrays, 
Haverah Park (UK)[l!l, Yakutsk (Russia) ^, Akeno and AGASA (Japan) Fly's Eye and HiRes (USA) 

we are probably very close to understanding the origin of UHECR. The forthcoming data of Auger detector (see Q 
for the first results) will undoubtedly shed more light on this problem. 

On the theoretical side we have an important clue to understanding the UHECR origin: the interaction of extra- 
galactic protons, nuclei and photons with CMB, which leaves the imprint on UHE proton spectrum, most notably in 
the form of the Greisen-Zatsepin-Kuzmin (GZK) [1, IH cutoff for the protons. 
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We shortly summarize the basic experimental results and the results of the data analysis, important for under- 
standing of UHECR origin (for a review see [13 )• 

(i) The spectra of UHECR are measured with good accuracy at 1 x 10^^ — 1 x 10^° eV, and these data have a power 
to reject or confirm some models. The discrepancy between the AGASA and HiRcs data at i? > 1 x 10^° eV might 
have the statistical origin (see jTlj and discussion in Section IIVC|) . and the GZK cutoff may exist. 

(ii) The mass composition at i? > 1 x 10^* eV (as well as below) is not well known (for a review see jT2l|). 
Different methods result in different mass composition, and the same methods disagree in different experiments. In 

principle, the most reliable method of measuring the mass composition is given by elongation rate (energy dependence 
of maximum depth of a shower, Xmax) measured by the fluorescent method. The data of Fly's Eye in 1994 favored 
iron nuclei at ~ 1 x 10"'^^ eV with a gradual transition to the protons at 1 x 10^^ eV. The further development of this 
method by the HiRes detector, which is the extension of Fly's Eye, shows the transition to the proton composition 
already at 1 x 10^^ eV [H [Ij. At present the data of HiRes [ij, [ij], HiRes-MIA flB'] and Yakutsk [16] favor the 
proton-dominant composition at ii^ > 1 x 10^* eV, data of Haverah Park |JJ do not contradict such composition at 
£: > (1 - 2) X 10^8 eV, while data of Fly's Eye [| and Akeno [13 agree with mixed composition dominated by iron. 
(in) The arrival directions of particles with energy E > A x 10^^ eV show the small-angle clusterin g w ithin the 
angular resolution of detectors. AGASA found 3 doublets and one triplet among 47 detected particles [l3| (see the 
discussion in [l^). In the combined data of several arrays [131 there were found 8 doublets and 2 triplets in 92 events. 
The stereo HiRes data ^2)] do not show small-angle clustering for 27 events at > 4 x 10^^ eV, maybe due to limited 
statistics. 

Small-angle clustering is most naturally explained in the case of rectilinear propagation as a random arrival of two 
(three) particles from a single source [l^]. This effect has been calculated in Refs. [1^ [13, [1^, [l^, [13, [H, [1^. In 
the last five works the calculations have been performed by the Monte Carlo (MC) simulations and results agree 
well. According to [11] the density of the sources, needed to explain the observed number of doublets is Ug — 
(1 — 3) X 10~^ Mpc"''. In [131 the best fit is given by ris ^ 1 x 10~^ Mpc~^ and the large uncertainties (in particular 
due to ones in observational data) are emphasized. 

(iv) Recently there have been found the statistically significant correlations between directions of particles with 
energies (4 — 8) x 10^^ eV and directions to AGN of the special type - BL Lacs [lO] (see also the criticism [sij and 
the reply [3^). 

The items (Hi) and (iv) favor rectilinear propagation of primaries from the point-like extragalactic sources, presum- 
ably AGN. However, the propagation in magnetic fields also exhibits clustering [25ll33l. IH] . 

The quasi-rectilinear propagation of ultra-high energy protons is found possible in MHD simulations [1^ of magnetic 
fields in large-scale structures of the universe (see however the simulations in [36l. [37j with quite different results). 

There are many unresolved problems in the field of Ultra High Energy Cosmic Rays, such as nature of primaries 
(protons? nuclei? or the other particles?), transition from galactic to extragalactic cosmic rays, sources and accelera- 
tion, but most intriguing problem remains existence of superGZK particles with energies higher than E ^ Ix 10^" eV. 
"The AGASA excess", namely 11 events with energy higher than 1 x 10^" eV, is still difficult to explain, though there 
are indications that it may have the statistical origin combined with systematic errors in energy determination (see 
section HVC|) . The AGASA excess, if it is real, should be described by another component of UHECR, most probably 
connected with the new physics: superheavy dark matter, new signal carriers, like e.g. light stable hadron and strongly 
interacting neutrino, the Lorentz invariance violation etc. 

The problem with superGZK particles is seen in other detectors, too. Apart from the AGASA events, there are five 
others: the golden Fly's Eye event with £' w 3 x 10^° eV, one HiRes event with E ^ 1.8 x 10^" eV and three Yakutsk 
events with E ^ 1 x 10^° eV. No sources are observed in the direction of these particles at the distance of order of 
attenuation length. The most severe problem is for the golden Fly's Eye event: with attenuation length ^att — 21 Mpc 
and the homogeneous magnetic field 1 nG on this scale, the deflection of particle is only 3.7°. Within this angle there 
are no remarkable sources at distance ~ 20 Mpc 3^] . 

In this paper we analyze the status of most conservative astrophysical solution to ultra-high energy cosmic ray 
problem, assuming that primary particles are protons or nuclei accelerated in extragalactic sources. In the first part 
of the paper (Sections II - V) we analyze the signatures of ultra-high energy protons propagating through CMB. We 
found that the dip, a spectral feature in energy range 1 x 10^^ — 4 x 10^^ eV, is well confirmed by observational data 
of AGASA, HiRes, Yakutsk and Fly's Eye detectors. In Sections VI -VII we discuss in the model-dependent way the 
transition from galactic to extragalactic cosmic rays and extragalactic sources: AGN and GRBs. 

We use in formulae for the flux throughout the paper, e.g. in equations (fTO)) . (fTTj) . p3|) . p4)) . (jSO]) and in Appendix 
lEl energy E measured in GeV, luminosity Lp - in GeV/s, emissivity (comoving energy density production rate) C - 
in GeV cm~^ s~^ and ii^min = 1 GeV, if not otherwise indicated. 
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II. ENERGY LOSSES AND THE UNIVERSAL SPECTRUM OF UHE PROTONS 



In this Section we present our recent calculations of energy losses for UHE protons interacting with CMB, and 
calculate the spectrum of protons, assuming the homogeneous distribution of sources in the space and continuous 
energy losses. The spectrum is calculated, using the conservation of the number of protons, interacting with CMB. 
Formally it does not depend on the mode of proton propagation (e.g. rectilinear or diffusive), and we shall discuss 
when this approximation is valid. The proton spectrum calculated in this way we call universal. 



A. Energy losses 

We present here the accurate calculations of energy losses due to pair production, p + jcmb p + e'^ + , and 
pion production, p + jcmb ^ N + pious, where ^cmb is a microwave photon. 

The energy losses of UHE proton per unit time due to its interaction with low energy photons is given by 



(1) 



where F is the Lorentz factor of the proton, Sr is the energy of background photon in the reference system of the 
proton at rest, eth is the threshold of the considered reaction in the rest system of the proton, (j{er) is the cross- 
section, f{er) is the mean fraction of energy lost by the proton in one P7 collision in the laboratory system, (see 
Fig. [22] in Appendix E)) . e is the energy of the background photon in the laboratory system, and n-y(e) is the density 
of background photons. 

For CMB with temperature T Eq. ([T]) is simplified 



1 dE cT 



E dt 27r2F2 
Further on we shall use the notation 



dera{er)f{er)er[^-\n l-cxp(--^) |. (2) 

^th 



for energy losses on CMB at present cosmological epoch, z — and T = Tq. For the epochs with red-shift z one has: 

l3{E,z) = il + zfpo[il + z)E]., (4) 
biE,z) = (1 + zfbo[{l + z)E]. (5) 

Another important physical quantity needed for calculations of spectra is the derivative dbo{E)/dE, which can be 
calculated as 

^ - -ME> + ^f,.M^,.)ns.) f_ ^ (6) 

exp ^ 2rTo J 

As one can see from Fig. [U dbo{E) / dE is numerically very similar to f3o{E), and for approximate calculations one 
can use /3o{E) values for both functions. 

From Eqs. ([T]) and ^ one can see that apart from cross-section the mean fraction of energy lost by the proton in 
laboratory system in one collision, /(e^) (see Eq. (jA2[) ). is the basic quantity needed for calculations of energy losses. 
The threshold values of these quantities are well known: 



_ 2me . _ '"-TT /_N 

/pair — 2^ I ™ ' /pion „, i ™ ' V ' / 



where /pair and /pion are the threshold fractions for p + j ^ p + + and p + 7 — > -I- tt, respectively. 

Pair production loss has been previously discussed in many papers. All authors directly or indirectly followed the 
standard approach of Ref. [s^ where the first Born approximation of the Bethe-Heitler cross-section with proton 
mass nip ~f 00 was used. In contrast to Ref. [39| . we use here the first Born approximation approach of Ref. (40| 
accounting for the finite proton mass. This allows us to calculate the average fraction of energy lost by the proton 
in laboratory system by performing a fourfold integration over invariant mass of electron-positron pair, Mx, over 
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FIG. 1: (a) UHE proton energy losses E~^dE/dt at z = (present work: curve 1; Berezinsky and Grigorieva (1988) 
curve 2; Stanev et al. 2000 [43| : black squares). The line 'red-shift' {Hq = 72 km/s Mpc) gives adiabatic energy losses. 
Note two important energies _Eoqi = 2.37 x lO'^* eV, where adiabatic and pair- production energy losses become equal, and 
-Eeq2 = 6.05 X 10^^ eV, where pair-production and photopion production energy losses are equal. The latter is one of the 
characteristics of the GZK cutoff, (b) The derivative dbo{E) / dE , where bo{E) = ~dE/dt at present epoch z — 0, (solid curve) 
in comparison with Po{E) = —E~^dE/dt (dashed curve). 



an angle between incident and scattered proton, and polar and azimuthal angles of an electron in the c.m.s. of the 
e"'"e~-pair (see Appendix A for further details). Calculating photopion energy loss we follow methods developed in 
papers [H, 44|. Total cross-sections are taken according to Ref. [45^. At low c.m.s. energy, Ec, we consider the binary 



reactions p -I- 7 +p, p-l-7 — > 7r+-|-n (they include the resonance p -I- 7 A+). Differential cross-sections of binary 
processes at small energies are taken from (46l . IZtI. |48| . At Ec > 4.3 GeV we assume the scaling behavior of differential 
cross-sections, the latter being taken from Ref. |49l. Isol. Isij . In the intermediate energy range we interpolate between 
the angular distribution of these two regimes with the cross-section being the difference between total cross-section 
and cross-section of the all binary processes. Angular distributions for this part of cross-section vary from isotropic 
at threshold to those imposed by inclusive pion photoproduction data at high energies. The overall differential cross- 
sections coincide with low-energy binary description and high-energy scaling distributions and join smoothly these 
two regimes in the intermediate region. 

The results of our calculations are presented in Fig. [1] in terms of E~^dE/dt as function of energy (curve 1). Also 
plotted is the derivative dbo{E)/dE (Fig.[T|3). This quantity is needed for calculation of differential energy spectrum. 
In Fig. [T] we plot for comparison the energy losses as calculated by Berezinsky and Grigorieva 1988 [llj (dashed curve 
2). The difference in energy losses due to pion production is very small, not exceeding 5% in the energy region relevant 
for comparison with experimental data {E < 10^^ eV). The difference with energy losses due to pair production is 
larger and reaches maximal value 15%. The results of calculations by Stanev et al. [i^ are shown by black squares. 
These authors have performed the detailed calculations for both aforementioned processes, though their approach is 
different from ours, especially for photopion process. Our new energy losses are practically indistinguishable from 
Stanev et al. [42| for pair production and pion production at low energies, and differ by 15-20% for pion production 
at highest energies (see Fig.[l]). Stanev et al. claimed that energy losses due to pair production is underestimated by 
Berezinsky and Grigorieva |41i | by 30-40%. Comparison of data files of Stanev et al. and Berezinsky and Grigorieva 
(see also Fig. [T]) shows that this difference is significantly less. Most probably, Stanev et al. scanned inaccurately the 
data from the journal version of the paper (ilj . 
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B. Universal spectrum 



To calculate the spectrum one should first of all evolve the proton energy E from the time of observation t = to 
(or z = 0) to the cosmological epoch of generation t (or red-shift z), using the the adiabatic energy losses EH{z) and 
h{E,z) given by Eq. g]): 



- dE/dt = EH{z) + (1 + zfho [(1 + z)E] , 



(8) 



where H{z) = i?o\/^^m(l + z)^ + is the Hubble parameter at cosmological epoch z, with Hq = 72 km/s Mpc, 
= 0.27 and f^A = 0.73 

We calculate the spectrum from conservation of number of particles in the comoving volume (protons change their 
energy but do not disappear). For the number of UHE protons per unit comoving volume, np{E), one has: 



/to 
dt Qg„,{Eg,t) dEg, 



(9) 



where t is an age of the universe, Eg = Eg{E,t) is a generation energy at age t, calculated according to Eq. ^ 
and QgcniEg,t) is the generation rate per unit comoving volume, which can be expressed through emissivity Co, the 
energy release per unit time and unit of comoving volume, at t = t^, as 



Qgc„(Sg, i) = £0(1 + z)"X<Zgo„(-Bg), 



(10) 



where (1 + z)"'' 
qgen{Eg) = Eg 



describes the possible cosmological evolution of the sources. In the case of the power-law generation, 
with normalization constant K = ^g — 2 for 7^ > 2 and K — (ln£',nax/-E'inin)^^ for 7^ = 2. We 
recall that in these formulae and everywhere below energies E are given in GeV, emissivity L in GeV cm^'^s^^ and 
source luminosity L in GeV s~^. 

From Eq. ^ one obtains the diffuse flux as 



Jp{E) ^ — CoKx 



An 







dt 


r 


dz 




'0 


dz 



dE„ 



(l + z)"Vn(i?,)^, 



(11) 



where \dt/dz\ = H^^{z)/{1 + z) and analytic expression for dEg/dE is given by Eq. (jB6p in Appendix [B] 

The spectrum (llip is referred to as universal spectrum. Formally it is derived from conservation of number of 
particles and does not depend on propagation mode (see Eq. But in fact, the homogeneity of the particles, 

tacitly assumed in this derivation, implies the homogeneity of the sources, and thus the condition of validity of 
universal spectrum is a small separation of sources. The homogeneous distribution of particles in case of homogeneous 
distribution of sources and inhomogeneous magnetic fields follows from the Liouville theorem (see Ref. [53|V 

Several effects could in principle modify the shape of universal spectrum. They include propagation in magnetic 
fields, discreteness in distribution of the sources, large-scale inhomogeneous distribution of sources, local source 
overdensity or deficit, and fluctuations in interaction. These effects will be studied in the next sections. With 
exception of energies beyond GZK cutoff and energies below 1 x 10^^ eV, the universal spectrum is only weakly 
modified by aforementioned effects. Here we shortly comment on role of magnetic fields. As numerical simulations 
(see e.g. [SJ, H^) show, the propagation of UHE protons in strong magnetic fields changes the energy spectrum (for 
physical explanation of this effect see [53|). The influence of magnetic field on spectrum depends on the separation 
of the sources, d. For uniform distribution of sources with separation d much less than characteristic lengths of 
propagation, such as attenuation length l^tt and the diffusion length Zdiff, the diffuse spectrum of UHECR is the 
universal one independent of mode of propagation (ssj . This statement has a status of the theorem. For the wide 
range of magnetic fields 0.1 - 10 nG and distances between sources d < 50 Mpc the spectrum at E > 10^^ eV is close 
to the universal one I56l. 



Modification factor 



The analysis of spectra is very convenient to perform in terms of the modification factor. 
Modification factor is defined as a ratio of the spectrum Jp{E) (see Eq. ([TT|)) with all energy losses taken into account, 
to the unmodified spectrum J™™(i?), where only adiabatic energy losses (red shift) are included. 



(12) 
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FIG. 2: Modification factor for the power-law generation spectra with 7g in the range 2.0 - 2.7. Curve t; = 1 corresponds to 
adiabatic energy losses, curves rjee - to adiabatic and pair production energy losses and curves r^tot - to all energy losses. The 
dip, seen at 1 x 10^* < E < ix 10^^ eV, has two flattenings: at low energy E ^ Ix 10^* eV and at high energy iJ fa 1 x 10^^ eV. 
The second flattening explains well the observed spectrum feature, known as ankle. 



For the power-law generation spectrum for 7^ > 2 without evolution one has 



c 





dt 


/ dz 




'0 


dz 



(13) 



The modification factor is a less model-dependent quantity than the spectrum. In particular, it should depend weakly 
on 7g, because both numerator and denominator in Eq. (|12p include E^'^" . In the next Section we consider the 
non-evolutionary case m — (see Section llV El for discussion of evolution). In Fig. [2] the modification factor is shown 
as a function of energy for two spectrum indices 7^ = 2.0 and jg = 2.7. As expected above, they do not differ much 
from each other. Note that by definition r]{E) < 1. 



III. SIGNATURES OF UHE PROTONS INTERACTING WITH CMB 



The extragalactic protons propagating through CMB have signatures in the form of three spectrum features: GZK 
cutoff, dip and bump. The dip is produced due to e"''e~-production and bump - by pile-up protons accumulated near 
beginning of the GZK cutoff. We add here the fourth signature: the second dip. 

The analysis of these features, especially dip and bump, is convenient to perform in terms of modification factor 
(4ll. [STj. For the GZK cutoff we shall use the traditional spectra. 



A. GZK cutoff 



The GZK cutoff [8|, |9| is most remarkable phenomenon, which describes the sharp steepening of the spectrum 
due to pion production. The GZK cutoff is a model-dependent feature of the spectrum, e.g. the GZK cutoff for a 
single source depends on the distance to the source. A common convention is that the GZK cutoff is defined for 
diffuse flux from the sources uniformly distributed over the universe. In this case one can give two definitions of 
the GZK-cutoff position. In the first one it is determined as the energy, Egzk ~ 4 x 10^^ eV , where the steep 
increase in the energy losses starts (see Fig. [T|). The GZK cutoff starts at this energy. The corresponding path length 
of a proton is Rgzk ~ {E~^dE /cdt)^^ w 1.3 x lO'^ Mpc. The advantage of this definition of the cutoff energy is 
independence of a spectrum index, but this energy is too low to judge about presence or absence of the cutoff in 
the measured spectrum. More practical definition is E1/2, where the fiux with cutoff becomes lower by factor 2 than 
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FIG. 3; -Ei/2 as characteristic of the GZK cutoff. The calculated integral spectra are multiplied to factor where 7 — 1 

found to fit the spectra in the interval 1 x 10^** — 3 x 10^^ eV. -B1/2 is shown by the vertical line. This value is valid for generation 
spectra with 2.1 < 7^ < 2.7. 



power-law extrapolation. This definition is convenient to use for the integral spectrum, which is better approximated 
by power-law function, than the differential one. In Fig. [3]the function E^'^^^\j{> E), where J(> E), the calculated 




E, eV E, eV 

FIG. 4: Predicted -E1/2 value in comparison with integral spectrum of Yakutsk array (left panel) and HiRes (right panel). One 
can see the agreement of the Yakutsk data with the theoretical value E1/2 ~ 5.3 x 10^^ eV (vertical line), while in the case of 
HiRes data this value is about 7 x 10^^ eV with large uncertainties due to difference between HiRes I and HiRes II data. HiRes 
collaboration found [58] £1/2 = (5.9lo:8) x lO" eV. 

integral diffuse spectrum, is plotted as function of energy. Note, that 7 > 7^ is an effective index of the power-law 
approximation of the spectrum modified by energy losses. For wide range of generation indices 2.1 < 7^ < 2.7 the 
cutoff energy is the same, E1/2 « 5.3 x 10^^ eV. The corresponding proton path length is R1/2 ~ 800 Mpc. In panel a) 
of Fig. d] El 12 is found from the integral spectrum of the Yakutsk array in the reasonable agreement with theoretical 
prediction. The HiRes data are shown in panel b). These data have large uncertainties, which prevent the accurate 
determination of i?i/2. However, they agree with the predicted value i?i/2 ~ 5.3 x 10^^ eV. In the recent paper [5l| 
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FIG. 5: UHECR differential spectra calculated with 73 = 2.7 and with different Emax are compared with AGASA data 
(left panel) and HiRes data (right panel). The curves (1),(2),(3) correspond to maximum acceleration energies Emax = 
3 X 10^° eV, 1 X 10^^ eV and 00, respectively. 





FIG. 6: Effect of local source overdensity and deficit on the shape of the GZK cutoff. The curve 1 is for the universal spectrum, 
the curves 2, 3 are for overdensity in the region of 30 Mpc with overdensity factor n/no equal to 2 and 3, respectively. The 
curves 5 and 6 show the deficit n/no = in the region of 10 Mpc and 30 Mpc, respectively. The calculated spectra are 
compared with the AGASA data (left panel) and with HiRes (right panel). The description of the AGASA excess, curve 4, 
needs unrealistic overdensity n/no = 10. 



the Hires collaboration found -E1/2 — (5.9j^Q g) x 10^^ eV in better agreement with the predicted value. 

In Fig. [5] the calculated universal spectra with the GZK cutoff are compared with AGASA and HiRes data 
While the HiRes data agree with the predicted GZK cutoff, the AGASA data show significant excess over this 
prediction at i? > 1 x 10^° eV. 

The GZK cutoff as calculated above has many uncertainties. 

The energy shape of the GZK feature is model dependent. The local excess of sources makes it flatter, and the deficit 
- steeper. The shape is affected by E^^^ (see Fig. [5]) and by fluctuations of source luminosities and distances between 
the sources. The cutoff, if discovered, can be produced as the acceleration cutoff (steepening below the maximum 
energy of acceleration in the generation spectrum). Since the shape of both, the GZK cutoff and acceleration cutoff, 
is model-dependent, it will be difficult to argue in favor of any of them, in case a cutoff is discovered. 
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To illustrate the effect of local overdensity/deficit of the UHECR sources we calculate here the UHECR spectra with 
different local ratios n/riQ, where n is the local density of the UHECR sources and rig is mean extragalactic source 
density. We use the various sizes of overdensity/deficit regions Roverd, equal to 10, 20 and 30 Mpc. The results of our 
calculations for the overdensity case are presented in Fig. [S] for jg = 2.7, to = and for four values of overdensity 
n/no equal to 1, 2, 3 and 10, assuming the size of overdensity region 30 Mpc (the results for Roverd = 20 Mpc are 
not much different). The spectra for the case of the deficit are shown by curves 5 and 6 in the both panels. The 
theoretical spectra shown in Fig. [S] illustrate uncertainties in the prediction of the shape of the GZK feature. 

An interesting question is whether the local overdensity of the sources can explain the AGASA excess. This problem 
has been already addressed in js^l for the realistic distribution of visible galaxies, considering them as the UHECR 
sources. Our calculations show (see left panel of Fig. (5]) that unrealistic overdensity n/np > 10 is needed to explain 
the AGASA excess. 

Both effects, Eyna,x and local source overdensity (deficit) affect weakly the shape of the GZK cutoff a.t E < 1 x 10^° eV. 
Thus, the precise measurements of the spectrum in this energy region, as well as measurement of E1/2, give the best 
test of the GZK cutoff. At higher energies theoretical predictions have large model-dependent uncertainties. 



B. Bump in the diffuse spectrum 



Protons do not disappear in the photopion interactions, they only loose energy and are accumulated near beginning 
of the GZK cutoff in the form of a bump. 

We see no indication of the bump in the modification-factor energy dependence in Fig. [2] As explained above, it 
should have been located at merging of riee{E) and r]tot{E) curves. The absence of the bump in the diffuse spectrum 
can be easily understood. The bumps are clearly seen in the spectra of the single remote sources (see left panel in 
Fig. [7]). These bumps, located at different energies, produce a flat feature, when they are summed up in the diffuse 
spectrum. This effect is illustrated by Fig. [7] (right panel). The diffuse flux there is calculated in the model where 
sources are distributed uniformly in the sphere of radius Rmax (or ^max)- When Zmax are small (between 0.01 and 0.1) 
the bumps are seen in the diffuse spectra. When radius of the sphere becomes larger, the bumps merge producing the 
flat feature in the spectrum. If the diffuse spectrum is plotted as E^Jp{E) this flat feature looks like a pseudo-bump. 



10-^ 



10^ 



§ 10 
« 10-^ 

UJ 

10' 



10' 
10'^ 



r=10Mpc 


30 


50 






100 






200 






500 . 






i' 1000 






1 1 1. . 




1 ■ 



10' 



10" 



10"* lO'^' 
E, eV 



10^ 




FIG. 7: Bumps for single sources at different distances (left panel) and disappearance of bumps in diffuse spectra (right panel). 
The calculations are performed for cosmological parameters j52j Hq — 72 km/sMpc, fim — 0.27 and f^A = 0.73, and numerically 
they are different from Ref. [4l|l. For the diffuse spectra (right panel) the sources are distributed uniformly in the sphere of 
radius 7?max, corresponding to ^max- The solid and dashed curves are for 73 — 2.7 and 7g = 2.0, respectively. The curves 
between Zmax ~ 0.2 and Zmax ~ 2.0 have Zmax = 0.3, 0.5, 1.0. 
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C. The dip 

The dip is more reliable signature of interaction of protons with CMB than GZK feature. The shape of the GZK 
feature is strongly model-dependent (see Section IIII Ap , while the shape of the dip is fixed and has a specific form 
which is difficult to mimic by other mechanisms, unless they have many free parameters. The protons in the dip 
are collected from the large volume with the radius about 1000 Mpc and therefore the assumption of the uniform 
distribution of sources within this volume is well justified, in contrast to the GZK cutoff, which strongly depends on 
local overdensity or deficit of the sources. The GZK cutoff can be mimicked by acceleration cutoff, and since the 
shape of the GZK cutoff is not rehably predicted, these two cases could be difficult to distinguish. 

The problem of identification of the dip depends on the accuracy of observational data, which should confirm the 
specific (and well predicted) shape of this feature. Do the present data have the needed accuracy? 

The comparison of the calculated modification factor with that obtained from the Akeno-AGASA data, using 
7g = 2.7, is given in Fig. [51 It shows the excellent agreement between predicted and observed modification factors for 
the dip. In Fig. [8] one observes that at E' < 1 x 10^* eV the agreement between calculated and observed modification 
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FIG. 8: Predicted dip in comparison with AGASA, HiRes, Yakutsk and Auger^T'] data. 

factors becomes worse and the observational modification factor becomes larger than 1. Since by definition ri{E) < 1, it 
signals the appearance of another component of cosmic rays, which is almost undoubtedly given by galactic cosmic rays. 
The condition 77 > 1 implies the dominance of the new (galactic) component, the transition occurs at i? < 1 x 10^* eV. 
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To calculate for the confirmation of the dip by Akeno-AGASA data, we choose the energy interval between 
1 X 10^^ eV and 4 x 10^^ eV (the energy of intersection of ijeeiE) and rjtatiE))- In calculations we used the Gaussian 
statistics for low-energy bins, and the Poisson statistics for the high energy bins of AGASA. It results in — 19.06. 
The number of Akeno-AGASA bins is 19. We use in calculations two free parameters: jg and the total normalization 
of spectrum. In effect, the confirmation of the dip is characterized by = 19.06 for d.o.f.=17, or x^/d.o.f.=1.12, 
very close to the ideal value 1.0 for the Poisson statistics. 

In the right-upper panel of Fig. [5] the comparison of modification factor with the HiRes data is shown. The 
agreement is also very good: x^ = 19.5 for d.o.f. = 19 for the Poisson statistics. The Yakutsk and Fly's eye data (not 
shown here) agree with dip as well. The Auger spectrum Q at this preliminary stage does not contradict the dip. 

The good agreement of the shape of the dip rj^e (E) with observations is a strong evidence for extra gala ctic protons 
interacting with CMB. This evidence is confirmed by the HiRes data on the mass composition Il4l|. While the 
data of the Yakutsk array [l^ and HiRes-MIA [l^ support this mass composition, and Haverah Park data [l| do 
not contradict it at £^ > (1 — 2) x 10^^ eV, the data of Akeno iJJl] and Fly's Eye [5] favor the mixed composition 
dominated by heavy nuclei. 

The observation of the dip should be considered as independent evidence in favor of proton-dominated primary 
composition in the energy range 1 x 10^^ — 4 x 10^^ eV. 



D. The second dip 



The second dip in the spectrum of extragalactic UHE protons appears at energy E = 6.3 x 10^^ due to interplay 
between pair production and photopion production. It is the direct consequence of energy £'cq2 = 6.05 x 10^^ eV, 
where energy losses due to e'^'e^-production and pion-production become equal (see Fig. [1]). This spectrum feature 
is explained as follows. 

The pion-production energy loss increases with energy very fast, and at energy slightly below i?cq2 e+e^-production 
dominates and spectrum can be with high accuracy described in continuous energy-loss approximation. At energy 
slightly higher than i?eq2 the pion-production dominates and the precise calculation of spectrum should be performed 
in the kinetic-equation approach. In this method the evolution of number of particles in interval dE is given by two 
compensating terms, describing the particle exit and regeneration due to collisions. The small continuous energy 
losses affect only the exit term and break this compensation, diminishing the flux. The exact calculations are given in 
Appendix [Dl The second dip is very narrow and its amplitude at maximum reaches ^ 10% (see Fig. lM)) . This feature 
can be observed by detectors with very good energy resolution, and it gives the precise mark for energy calibration 
of a detector. It can be observed only marginally by the Auger detector. 



IV. ROBUSTNESS OF THE DIP PREDICTION 



We calculated the dip for the universal spectrum, i.e. for the case when distances between sources are small enough, 
and the spectrum does not depend on the propagation mode. In this Section we shall study stability of the dip relative 
to other possible phenomena, namely, discreteness in the source distribution, propagation in magnetic fields etc. We 
shall consider also some phenomena related to existence of the dip. 



A. Discreteness in the source distribution 



As it follows from analysis of the small-scale anisotropy (see [12, Hj, [Ij, [H, Hal), the average distance between 
UHECR sources is d ^ 30—50 Mpc [27l. [28] . Such discreteness affects the spectrum, especially at highest energies, when 
energy attenuation length is comparable with d. In this subsection we demonstrate the stability of the dip relative 
to discreteness of the sources. We illustrate the effect of discreteness by an example of UHE protons propagating 
rcctilinearly from sources located in the vertices of a 3D cubic lattice with spacing d. Positions of sources are given 
by coordinates x — id, y — jd and z ~ kd, where k — 0, ±1, ±2 . . .. The observer is assumed a,tx = y = z = 
with no source there. The diffuse flux for the power-law generation spectrum oc Eg is given by summation over all 
vertices. The maximum distance is defined by the maximum red-shift Zmax- Then the observed flux is given by 

(7g - ^)^od Eg{E, z.jkY^o dEg 

' (47r)2 + fc2)(i + ^^^.^) dE ' ^'^^ 
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FIG. 9: Proton spectra for rectilinear propagation from discrete sources. Sources are located in vertices of 3D cubic grid with 
spacing d = 60, 40, 20, 10, 5 and 1 Mpc. The calculations are performed for Zmax — 4, Emax = 1 x 10^^ eV and 79 — 2.7. 

where Co = Lp/(P is emissivity, z^jfe is the red-shift for a source with coordinates i,j,k, and factor (1 + Zijk) takes 
into account the time dilation. 

The calculated spectra for d — 1,5, 10, 20, 40 and 60 Mpc are shown in Fig.[n]in comparison with the AGASA-Akeno 
data. In calculations we used E'max = 1 x 10^^ eV, m = (no evolution), z^ax = 4 and 7^ = 2.7. Emissivity Cq is 
chosen to fit the AGASA data. One can see that discreteness in the source distribution alTects weakly the dip, but 
the effect is more noticeable for the shape of the GZK cutoff. 

With d decreasing, the calculated spectra regularly converge to the universal one, as it should be according to 
propagation theorem [s^ . This theorem ensures also that the spectra from Fig. [14] are valid for the case of weak 
magnetic field, when the diffusion length Zdiff > d. 

B. Dip in the case of diffusive propagation 

The dip, seen in the universal spectrum, is also present in the case of diffusive propagation in magnetic field f56| . 
The calculations are performed for diffusion in random magnetic fields with the coherent magnetic field Bq = 1 nG 
and up to Bq = 100 nG on the basic scale Ic = 1 Mpc. The calculated spectrum is shown in Fig. [10] for the case 
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FIG. 10: Diffusive energy spectrum in the case of Bo = 1 nG, Ic = 1 Mpc and for the diffusion regimes: Kolmogorov (continuous 
line), Bohm (dashed line) and D{E) oc (dotted line). The separation between sources is rf = 50 Mpc and the injection 
spectrum index is 7^ = 2.7. The universal spectrum (dash-dotted line) and the AGASA-Akeno data are also shown. 

{Bq,Ic) = (1 nG, 1 Mpc) and distance between sources d = 50 Mpc. The critical energy, where diffusion changes its 
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regime is i5b « 1 x 10"'^^ eV. The spectra are shown for three different regimes a.t E < Eb as indicated in Fig. [TOl The 
universal spectrum is also presented. One can see that the dip agrees well with universal spectrum and observational 
data, while the shape of the GZK cutoff differs considerably from the universal spectrum. 

These calculations demonstrate stability of the dip relative to changing of the propagation mode, and sensitivity of 
the GZK cutoff to the way of propagation. 



C. Energy calibration of the detectors using the dip 



Since the position and shape of the dip is robustly fixed by proton interaction with CMB, it can be used as energy 
calibrator for the detectors. We use the following procedure for the calibration. Assuming the energy-independent 
systematic error, we shift the energies in each given experiment by factor A to reach minimum for comparison 
with the calculated dip. The systematic errors in energy determination of existing detectors exceed 20%, and it 
determines the expected value of A. The described procedure results in A^ = 0.9, Aya = 0.75 and Xm = 1.2 for 
AGASA, Yakutsk and HiRes detectors, respectively. After this energy calibration the fluxes in all experiments agree 
with each other. First we consider the AGASA and HiRes data. There are two discrepancies between these data (see 
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FIG. 11: Spectra of Akeno-AGASA, HiRes and Yakutsk before and after energy calibration by the dip. The spectra with 
energy shift are shown in the right panels. The energy shifts needed for the best fit of the dip are AAg ~ 0.9 , Ahi ~ 1.2 and 
Ava = 0.75 for AGASA, HiRes and Yakutsk, respectively. 



the upper-left panel of Fig. [TT|) : one is described by factor 1.5 - 2.0 in energy region 1 x 10^^ — 8 x 10^^ eV, and the 
second - at > 1 x 10^° eV. In Fig. [TT] the spectra of Akeno-AGASA and HiRes are shown before and after the 
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energy calibration. One can see the good agreement of the cahbrated data at i5 < 1 x 10^ eV and their consistency 
at i? > 1 X 10^*^ eV. This result should be considered together with calculations [gOl, where it was demonstrated 
that 11 superGZK AG AS A events can be simulated by the spectrum with GZK cutoff in case of 30% error in energy 
determination. We may tentatively conclude that existing discrepancy between AGASA and HiRes spectra at all 
energies are due to systematic energy errors and statistics. 

In Fig.[TT]the energy spectra are shown for AGASA and Yakutsk spectra before and after energy calibration. Again, 
the best fit to the dip shape results in excellent agreement in the absolute values of fluxes. 

The agreement between spectra of all three detectors after energy calibration by the dip confirms the dip as the 
spectrum feature produced by interaction of the protons with CMB, and demonstrates compatibility of fluxes measured 
by AGASA, HiRes and Yakutsk detectors. 



D. Dip and extragalactic UHE nuclei 

The proton dip has very good agreement with observations (see Fig. [8]). However, in all astrophysical sources the 
nuclei must be also accelerated to the energies, naively, Z times higher than that for protons. Do UHE nuclei in 
primary flux upset the good agreement seen in Fig. [8]/ This problem has been recently considered in [6l|, HI, HI] (for 
study of propagation UHE nuclei through CMB see [s^, [fi^ EHl)- The presence of nuclei in primary extragalactic 
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FIG. 12: Modification factors for helium and iron nuclei in comparison with that for protons. Proton modification factors are 
given by curves 1 and 2. Nuclei modification factors are given by curves 3 (adiabatic and pair production energy losses) and 
by curves 4 (with photodisintegration included). 



spectrum modifles the dip [61| . In Fig. [12] the dips for helium and iron nuclei are shown in comparison with that 
for protons. From this flgure one can see that presence of 15 - 20 % of nuclei in the primary flux breaks the good 
agreement of proton deep with observations. The modiflcation factor for cosmic rays composed of protons and nuclei 
with the fraction — Qa{E)/Qp{E), where Q{E) is generation function for nuclei (A) and protons (p), can be 
easily calculated as 



r?tot(£;) 



1 + A 



(15) 



The fraction A is a model-dependent value, which depends on ratio of number densities of gas components ua/tih 
in media, where acceleration operates, on ionization of the gas and on the injection mechanism of acceleration [67j . 
Besides, the UHE nuclei can be destroyed inside the source or in its vicinity [63]. 

The strongest distortion of proton modification factor is given by helium nuclei, for which nRc/nu ~ 0.08 cor- 
responding to the helium mass fraction Yp = 0.24. In Fig. [12] the modification factors are shown for the mixed 
composition of protons and helium with the mixing parameter A = 0.1 (left panel) and A = 0.2 (right panel). One 
can judge from these graphs about allowed values of mixing parameter A. If agreement of the proton dip with obser- 
vations is not incidental (the probability of this is small according to small x^/d.o.f.). Fig. [T^] should be interpreted 
as indication to possible acceleration mechanism 67]. 
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E. Dip and cosmological evolution of the sources 

The cosmological evolution of the sources, i.e. increase of the luminosities and/or space densities with red-shift 
z, is observed for many astronomical populations. The evolution is reliably observed for star formation rate in the 
normal galaxies, but this case is irrelevant for UHECR, because neither stars nor normal galaxies can be the UHECR 
sources due to low cosmic-ray luminosities Lp and maximum energy of acceleration -Emax- AGN, which satisfy 
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FIG. 14: The dip in evolutionary models in comparison with the AGASA data. The parameters of evolution used in the 
calculations for curves 1 and 2 are similar to those observed for AGN. The curve 3 is for m = Q. 

these requirements, also exhibit the evolution seen in radio, optical and X-ray observations. The X-ray radiation is 
probably most relevant tracer for evolution of UHECR because both radiations are feed by the energy release provided 
by accretion to a massive black hole: X-rays - through radiation of accretion disk, and UHECR ~ through acceleration 
in the jets. According to recent detailed analysis in [68| and [69| the evolution of AGN seen in X ray radiation can 
be described by factor (1 -I- z)™ up to Zc ~ 1.2 and is saturated at larger z. In [6l| the pure luminosity evolution 
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and pure density evolution are allowed with m = 2.7 and m = 4.2, respectively, and with Zc « 1.2 for both cases. 
In [69| the pure luminosity evolution is considered as preferable with m = 3.2 and Zc = 1-2. These authors do not 
distinguish between different morphological types of AGN. It is possible that some AGN undergo weak cosmological 
evolution, or no evolution at all. BL Lacs [301 , which are important as potential UHECR sources, show no signs of 
positive cosmological evolution [70] . 

In case of UHECR there is no need to distinguish between luminosity and density evolution, because the diffuse 
flux is determined by the emissivity, which includes both luminosity and space evolution, as it follows from Eq. (jlOp . 

In Fig. [T3] we present the calculated dip spectrum in evolutionary models, inspired by the data cited above. For 
comparison we show also the case of absence of evolution m = 0, as can be valid for BL Lacs. From Fig. [14] one 
can see that the spectra with evolution up to > 1 can explain the observational data down to (a few) x 10^^ eV 
and even below, in accordance with early calculations [ti], [t^ (see [tI] for recent analysis). However, for any 
reasonable magnetic fields protons with these energies have small diffusion lengths and the spectrum acquires the 
diffusion 'cutoff at energy Eb = 1 x 10^^ eV (see Section HVBl) . 

We conclude that for many reasonable evolution regimes the dip agrees with observational data as well as the 
non-evolutionary case m = 0. 



V. ROLE OF INTERACTION FLUCTUATIONS 



UHE proton spectrum is affected by fluctuations in the photopion production. These fluctuations may change the 
proton spectrum only at energy substantially higher than E — Ax 10^^ eV. At this energy the half of energy losses is 
caused by e~^e~ production which does not fluctuate. Up to energy 1 x 10^" eV the photoproduction of pions occurs 
at the threshold in collisions with photons from high-energy tail of the Planck distribution, and fraction of energy 
lost does not fluctuate, being flxed by the threshold value. Indeed, for = 1 x 10^° eV the minimal energy of CMB 
photon needed for pion production is e = 3 x 10~^ eV to be compared with energy of photon in the Planck distribution 
maximum = 3.7 x 10~^ eV. The only fluctuating value is the interaction length. 

The noticeable effect of fluctuations is expected for protons with energies E > 1 x 10^" eV. 

As it is well known (tsI . [76j . the kinetic equations give an adequate method to account for the fluctuations in 
interaction. Neglecting the conversion of proton to neutron (neutron decays back to proton with small energy loss) 



the kinetic equation for UHE protons with adiabatic energy losses and with p- 
scattering in collisions with CMB photons can be written down as follows: 



■7 ■ 



p + e'^ + e and p + ^ ^ N +pions 



On {E t] d 

' ' =-3H{t)np{E,t) + — {[Hit)E + bpa;riE,t)]np{E,t)} 



dt 



P{E,t)np{E,t) + 



(16) 



dE'P{E\ E, t)np{E\ t) + QgeniE, t), 



where np{E,t) is the number density of UHE protons per unit energy, QgeniE,t) is the generation rate, given by 
Eq. ([TO|l with m = and H{t) is the Hubble parameter. The first term in the r.h.s. of Eq. ([T6l) describes expansion 
of the universe. The energy loss hpair{E) due to e+e~-pair production is treated as continuous energy loss. The 
photopion collisions are described with help of probability P{E, t) for proton exit from energy interval (i?, E + dE) 
due to p7-collisions and with help of their regeneration in the same energy interval described by probability P{E' , E, t). 
These two probabilities describe fluctuations in the interaction length and in fraction of energy lost in the interaction; 
the interaction length is equivalent in this picture to time of proton exit from energy interval dE. The exit probability 
P{E,t) due to collisions with CMB photons with temperature T can be written as: 



P{E,t) 



cT 



Att^E^ 



dE,E,{El - ml)a{E,) In 



1 — exp 



ml 



AET 



(17) 



where E^. is the total c.m.s. energy of colliding proton and photon, E™^"^ — itIt^ + irip, a{Ec) is the photopion cross- 
section and T is the CMB temperature at cosmological epoch t. 

Similarly, in regeneration term of Eq. (|16p the probability for a proton with energy E' to produce a proton with 
energy E is given by 



P{E',E,t)^- 



cT 
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(18) 
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FIG. 15: Modification factor for the power-law generation spectra with 75 — 2.7 and Zma^ = 5. Upper dotted fine corresponds 
to the case of adiabatic energy loss. Curve 'e'^e~' gives the modification factor for the case of adiabatic and pair-production 
energy losses. Curve 'cont.loss' corresponds to total energy losses in continuous loss approximation. The 'fluct.' curve describes 
the case of numerical solution to the kinetic equation (|16p . The left and right panels correspond to iJmax = 1 x 10^'^ eV and 
-Emax = 1 X 10^^ eV, respectively. 



where E' and E are the energies of primary and secondary protons, respectively, and x = E/E' . The minimum value 
of the allowed c.m.s. energy in this case is given by 

Er''{x)^[ml/il-x) + ml/x]'/'. (19) 

This bound corresponds to the process with minimum invariant mass, namely to p + 7 — > tt*^ + p. 

We have solved Eq. numerically. The calculated spectrum Jp{E) ~ (c/47r)np(i?, ip), where tg is the age of 
the universe at red-shift z = 0, is presented in Fig. [13 as modification factor 7]{E) = Jp{E)/ Jp'^'^{E) for generation 
spectrum with 7^ = 2.7 and E^ax = 1 x 10^^ eV (left panel) and E^ax = 1 x 10^^ eV (right panel). For comparison 
the modification factors for universal spectrum with continuous energy losses is also shown. The difference in these 
two spectra at highest energies must be due to fluctuations in energy losses, though formally we have to say that 
this is the difference between solution to kinetic equation (jl6p and the continuous energy loss approximation. For 
Emax = 1 X 10^'^ eV one can see the difference in the spectra about 25% at highest energies and tiny difference above 
intersection of rj^e and rytot curves. For E^ax = 1 x 10^^ eV the difference is small. 

Note, that modification factors do not vanish at E-^^-^ even when generation function goes abruptly to zero, since 
both solutions vanish keeping the same value of ratios i]{E) = Jp{E) / J^nmiE) . It is easy to demonstrate analytically 
that ratio of flux in continuous loss approximation Jcont{E) to unmodified flux Junm{E), given by Eq. ([13]) tends to 
^^o//?coii(£'max), when E — > -Emax. From Fig. [T5| one can see that this ratio coincides exactly with our numerical 
calculations (e.g. the analytical value is 2.45 x 10^'' for i?niax = 1 x 10^'^ eV), and this gives a proof that our numerical 
calculations are correct. 

The effect of interaction fluctuations is usually taken into account with help of Monte-Carlo simulations. The 
method of kinetic equation corresponds to averaging over large number of Monte-Carlo simulations, and if all other 
assumptions are the same, the results must coincide exactly. These assumptions include i?max and parameters of 
P7-interaction. However, the existing Monte-Carlo simulations in most cases include some other assumptions in 
comparison with kinetic equations, which modify spectrum stronger than interaction fluctuations. One of them is 
discreteness in the source distribution (in kinetic equations the homogeneous distribution is assumed), the other is 
fluctuations of distances to the nearby sources. 

It is possible however to make the comparison with Monte-Carlo simulations for homogeneously distributed sources 
and using identical interaction model. Such a comparison is discussed in Appendix [Cl 

As to results presented here, it is necessary to emphasize that the difference of kinetic-equation solution and 
continuous-energy-loss approximation presented in Fig. [T51 includes fluctuations, but not only fluctuations. The tran- 
sition of kinetic equation to continuous-energy-loss equation depends on some other conditions which can fail. What 
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the presented calculations demonstrate (see Fig. 1151) is that continuous-energy-loss approximation describes with very 
good accuracy more reliable kinetic-equation solution, which in particular includes interaction fluctuations. 



VI. TRANSITION FROM EXTRAGALACTIC TO GALACTIC COSMIC RAYS 

In the analysis above we obtained several indications that transition from extragalactic to galactic cosmic rays 
occurs at i? « 1 X 10^^ eV. These evidences are summarized in Fig. [TH 

The predicted spectrum above 1 x 10^^ eV describes perfectly well the observed spectra: see the modification factor 
in the upper-left panel of Fig. [TBI compared with AGASA data and in Fig. [8] with HiRes. However, at i? < 1 x 10^^ 
experimental modification factor becomes > 1, in contrast to definition rj < 1. It signals the appearance of a 
new component, which can be nothing but the galactic cosmic rays. In the right panel the spectrum for rectilinear 
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FIG. 16: Appearance of transition energy ~ 1 x 
10^* eV in the modification factor compared with 
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propagation from the sources with different separation d and with 7^ = 2.7 is compared with AGASA data. One can 
see that at i? < 1 x 10^^ eV the calculated extragalactic spectrum becomes less than that observed. 

In the lower panel the similar comparison is shown for diffusive propagation with different diffusion regimes at lower 
energies [5^ [73| . The random magnetic field with basic scale Ic = ^ Mpc and magnetic field on this scale Bq = 1 nG 
is assumed. The dash-dotted curve (universal spectrum) corresponds to the case when the separation between sources 

In all cases the transition from extragalactic to galactic component begins at w 1 x 10^^ eV, with index b for 
"beginning" . 

What is the reason of this universality? We study the transition, moving from high towards low energies. Eh is the 
beginning of transition (or its end, if one moves from low energies). Ei, is determined by energy E^qi ~ 2.37 x 10^^ eV, 
where adiabatic and pair-production energy losses become equal. The quantitative analysis of this connection is given 
in [5^. We shall give here the semi-quantitative explanation. 



19 



The flattening of the spectrum occurs at energies E < Ei,, where Ei, = Eeq/{1 + Zes)'^ and ZeS should be estimated 
as red-shift up to which the main contribution to unmodified spectrum occurs. The simphfied analytic estimate for 
7g = 2.6 — 2.8 gives 1 + z^b w 1.5 and hence i?h w 1 x 10^^ eV. In fact, the right and lower panels of Fig. [T6l present 
the exact calculations of this kind. 

In experimental data the transition is searched for as a feature started at some low energy E2kn - the second knee. 
Its determination depends on experimental procedure, and all we can predict is i?2kn < Eb- Determined in different 
experiments E2kn ^ (0.4 - 0.8) x 10^^ eV. 

The transition at the second knee appears also in the study of propagation of cosmic rays in the Galaxy (see e.g. 

Being thought of as purely galactic feature, the position of the second knee in our analysis appears as direct 
consequence of extragalactic proton energy losses. 




E, GeV E, GeV 

FIG. 17: Transition from extragalactic to galactic cosmic rays in the second knee (left panel) and ankle (right panel) models. 
In the left panel are shown: KASCADE total spectrum, which above the iron knee -Bpe is composed mostly by iron nuclei 
('gal.Fe' curve), below Et the extragalactic proton spectrum ('extr.p' curve) is calculated for diffusive propagation (see Section 
IVIIB|) and Etr is the energy of transition from galactic cosmic rays to extragalactic protons. The dot-dash line shows power-law 
extrapolation of low-energy KASCADE spectrum. In the right panel extragalactic proton spectrum is calculated for generation 
spectrum oc E~^, while galactic spectrum (curve 'gal.CR') is taken as difference between the observed total spectrum and 
calculated spectrum of extragalactic protons. 

The transition at the second knee is illustrated by Fig. [iTl The clue to understanding of this transition is given by 
the KASCADE data [8l|, [S^l • They confirm the rigidity model, according to which position of a knee for nuclei with 
charge Z is connected with the position of the proton knee Ep as Ez = ZEp. There are two versions of this model. 
One is the confinement-rigidity model (bending above the knee is due to insufficient confinement in galactic magnetic 
field), and the other is acceleration- rigidity model (-Emax is determined by rigidity). In both models the heaviest nuclei 
(iron) start to disappear at E' > Epc = 6.5 x 10^^ eV, if the proton knee is located at Ep « 2.5 x 10^^ eV. The shape 
of the spectrum above the iron knee (E > Efo) is model dependent, with two reliably predicted features: it must be 
steeper than the spectrum below the iron knee [E < E-pc), shown by the dash-dot curve, and iron nuclei must be 
the dominant component there (see Fig. I17p . The high energy part of the spectrum has a characteristic energy Et, 
below which the spectrum becomes more fiat, i.e. drops down when multiplied to E^'^ (see Fig. I17p . This part of the 
spectrum is shown for the diffusive propagation described in section I VII Bl These two falling parts of the spectrum 
inevitably intersect at some energy Etr, which can be defined as transition energy from galactic to extragalactic cosmic 
rays. The 'end' of galactic cosmic rays Epe = 6.5 x 10^^ eV and the beginning of full dominance of extragalactic 
component E^ ^ Ix 10^^ eV differ by an order of magnitude. Note, that power-law extrapolation of the total galactic 
spectrum, shown by dot-dash line, beyond the iron knee Epo bas no physical meaning in the rigidity models and must 
not be discussed. 

The second-knee transition gives an alternative possibility in comparison with ankle-transition hypothesis known 
from end of 1970s. It is inspired by ffattening of the spectrum at Ea ~ 1 x 10"'^^ eV seen in the AGASA and Yakutsk 
data (left panels in Fig. [5]) and possibly at (0.5 — 1) x 10"^^ eV in the Hires data (the right panel in Fig. [5]). Being 
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multiplied to factor E'^-^, as in Fig.[T71 the ankle transition looks very similar to that at the second knee. Note that 
in the latter case the ankle is just an intrinsic part of the dip. 

The ankle transition has been recently discussed in Refs. [H, [H, [H, [H, [s^l • 

In the ankle model it is assumed that galactic cosmic ray spectrum has a power-law shape oc E~'^ from the proton 
knee Ep w 2.5 x 10^'^ eV to about Ea ^ I x 10^^ eV where it becomes steeper and crosses the more flat extragalactic 
spectrum (see the right panel of Fig. [T7|) . The ankle transition in Fig. [IT] is shown for extragalactic proton spectrum 
with generation index jg = 2, while galactic spectrum, given by curve "gal.CR" is calculated as difference of the 
observed total spectrum and calculated extragalactic proton spectrum. 

The ankle model has the problems with galactic component of cosmic rays. The spectrum at 1 x 10^* — 1 x 10^^ eV 
is taken ad hoc to fit the observations, while in the second knee model this part of the spectrum is calculated 
with excellent agreement with the data. In the rigidity models the heaviest nuclei (iron) start to disappear at 
E> Ef^ = 6.5 x 10^^ eV. How the gap between 1 x lO^^ eV and 1 x lO^^ eV is filled? 

Galactic protons start to disappear at > 2.5 x 10^^ eV. Where they came from at _E > 1 x 10^^ eV to be seen 
e.g. in the Akeno detector with fraction ^ 10%? 

The ankle model needs acceleration by galactic sources up to 1 x 10^^ eV (at least for iron nuclei), which is difficult 
to afford. The second knee model ameliorates this requirement by one order of magnitude. 

The second-knee model predicts the spectrum shape down to 1 x 10^^ eV with extremely good accuracy (x^/d.o.f.= 
1.12 for Akeno-AGASA and x^/d.o.f.= 1.03 for HiRes). In the ankle model one has to consider this agreement as 
accidental, though such hypothesis has very low probability, determined by cited above. As an alternative the 
ankle model-builders can suggest only hopes for future development of galactic propagation models to be as precisely 
calculated as the dip. 



VII. ASTROPHYSICAL SOURCES OF UHECR 



In the sections above we have performed the model-independent analysis of spectra of extragalactic protons interact- 
ing with CMB. We have calculated the features of the proton spectrum assuming the power-law generation spectrum 
(X valid at > 1 X 10^^ eV, and compared predicted features with observations. We found that proton dip, a 

model-independent feature at energy between 1 x 10^^ eV and 4 x 10"'^^ eV, is well confirmed by observations. Only two 
free parameters are involved in fitting of observational data: -fg = 2.7 (the allowed range is 2.55 - 2.75) and the flux 
normalization constant. The various physical phenomena included in calculations, such as discreteness in the source 
distribution, the different modes of propagation (rectilinear and diffusive), cosmological evolution with parameters 
similar to AGN evolution, fluctuations in pj interaction etc, do not upset this agreement. 

The transition of extragalactic to galactic cosmic rays is also discussed basically in model-independent manner. 

In this Section we shall discuss the models: realistic energy spectra, the sources and the models for transition from 
extragalactic to galactic cosmic rays. 

The UHECR sources have to satisfy two conditions: they must be very pow erful and must accelerate particles to 
large i?max ^ 1 x 10^^ eV. There is one more restriction [g^, IH, [13, [25l. l26l. [27l . \2§^ . coming from observation of 
small-scale clustering: the space density of the sources should be (1 — 3) x 10~^ Mpc""^ probably with noticeable 
uncertainty in this value. Thus, these sources are more rare, than typical representatives of AGN, e.g. the Seyfert 
galaxies, whose space density is ~ 3 x 10""* Mpc~'^. The sources could be the rare types of AGN, and indeed the 
analysis of [13] show statistically significant correlation between directions of par ticles with energies (4 — 8) x 10^^ 
eV and directions to AGN of the particular type - BL Lacs (see also criticism f3l'| and reply [1^). The acceleration 
in AGN can provide the maximum energy of acceleration up to ~ 10^^ eV for non-relativistic shock acceleration (see 
e.g. Hi). 

The relativistic shock acceleration can occur in AGN jets. Acceleration to i^max ~ 10^^ eV in the AGN relativistic 
shocks is questionable (see discussion below). 

Gamma Ray Bursts (GRBs) are another potentially possible sources of UHECR [8^[93, [9l|. They have very large 
energy output and can accelerate particles up to ^ 10^^ eV [93,[9l|. These sources have, however, the problems with 
explaining small-angle anisotropy and with energetics (see discussion below). 



A. Spectra 

The assumption of the power-law generation spectrum with 7^ = 2.7 extrapolated to E'min ~ 1 GeV results in too 
large emissivity required for observed fiuxes of UHECR. To avoid this problem the broken generation spectrum has 
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been suggested in Refs. [tHIq^: 

^g^'^^^) { a at E>eI ' 

where Qgcn{E) is the generation function (rate of particle production per unit of comoving volume), defined by Eqs. ^ 
and (jlOp. and E'c was considered as a free parameter. 

Recently it was demonstrated that broken generation spectrum can naturally emerge under most reasonable physical 
assumption. In Ref. [93l | it has been argued that while spectrum E^^ is universal for non-relativistic shock acceleration, 
the maximum acceleration energy i?max is not, being dependent on the physical characteristics of a source, such as 
its size, regular magnetic field etc. Distribution of sources over i^max results in steepening of the generation function, 
so that the distribution of the sources dn/d-Emax oc E~^^ explains the observational data, if /? = 1.5 — 1.6 [93|. 

In Appendix[E]we address the generalized problem: what should be the distribution of spectral emissivity over i?max 
to provide the generation function with the broken spectrum. We use there the notation e = -Emax and introduce the 
spectral emissivity C{e) = ns{e)Lp(e), where Lp{e) is particle luminosity of a source and ns{e) = ns{e^ Lp(e)) is the 
space density of the sources. The total emissivity is given by £o = / C{e)de. Distribution of spectral emissivity C{e) 
over maximal energies e determines the energy steepening of the generation function Q^cn{E) at energy E — Emin in 
the distribution. This function is calculated analytically for arbitrary £(£), assuming that e is confined to the interval 
(emin , £max)- It is demonstrated that QgcniE) can be the power-law function cx E~''^ exactly, only if C(e) distribution 
is power-law, too (cx e^^). For the source generation function qgcn{E) cx g± E < e, the generation index 

in the interval emin < E < emax is found to be 7g = 1 + a + /S, including the case a = 0. The steepening of the 
generation spectrum from 2 -I- a to 7g occurs approximately at energy Ec = £min- At energy E > emax the spectrum 
is suppressed as exp(— i?/emax) or stronger. 

In the applications we are interested in two cases: a = (non-relativistic shocks) and a = 0.2 — 0.3 (ultra-relativistic 
shocks). In the latter case the term E^^ in Eq. (^0)1 should be substituted by The energy Ec — Smin in 

Eq. ([^D|) is considered as a free parameter. 

B. Active Galactic Nuclei 

The AGN as sources of UHECR meet the necessary requirements: (i) to accelerate particles to £^max ^ 10^"'^ eV, 
(ii) to provide the necessary energy output and (Hi) to have the space density ris ^ (1 — 3) x 10~^ Mpc^"^, required 
by small-scale clustering. We shall discuss below these problems in some details. 

1. Acceleration and spectra 

The flow of the gas in AGN jet can be terminated by the non-relativistic shock which accelerates protons or nuclei 
in the radio lobe up to i?max ~ 10^^ eV with spectrum cx E^^ 88]. 

In some cases the observed velocities in AGN jets are ultra-relativistic with Lorentz factor up to F ~ 5 — 10. It 
is natural to assume there the existence of internal and external ultra-relativistic shocks. Acceleration in relativistic 
shocks relevant for UHECR has been recently studied in Refs. [13, HH, [H, HE S IS H HI- The acceleration 
spectrum is cx E~'^!' and in the case of isotropic scattering of particles upstream and downstream, the spectrum index 
is 7g = 2.23 ± 0.01 |95|]. However, recently it was understood that this result depends on scattering properties of the 
medium [13, [9§], and the spectrum can be steeper. In the regime of large angle scattering in [o^l, 7g = 2.7 was found 
possible for the shock with velocity u ~ (0.8 — 0.9)c and compression ratio r = 2. In Monte-Carlo simulation [9^ it 
is demonstrated that effect of compression of upstream magnetic field results in increasing of 7g up to the limiting 
value 2.7 in ultra-relativistic case Fgh ^ 1- 

The maximal energy of acceleration i?max is a controversial issue. While in most works (very notabl y 1941 ) it is 
obtained that E'max cannot reach ^ eV for all realistic cases of relativistic shocks, the authors of Ref. [lOOf argue 
against this conclusion. 

We shall divide a problem with i?max into two: the reliably estimated energy gain in relativistic shocks and model- 
dependent absolute value of i?max- The energy gain in the first full cycle of particle reflection upstream-downstream- 
upstream (u — > d — > u) is about F^j^. The next reflections are much less effective, as was first observed in [93|: a 
particle lives a short time in the upstream region, before it is caught up by the shock. As a result, a particle is deflected 
by upstream magnetic field only to a small angle, and thus it occurs in the downstream region with approximately 
the same energy as in the first cycle. Then the energy upstream will be also almost the same as in the first cycle. 
According to Monte-Carlo simulation [9^ the average energy gain per each successive u ^ d ^ u cycle is only 1.7 (see 
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Fig. 4 in [931 with clear explanation). With these energy gains (~ F^j^ in the first u ^ d ^ u cycle and with '--^ 2 in 
each successive cycle) it is not possible to get i?max 10^^ e V in the conservative approach for AGN and GRBs, but 
there are some caveats in this conclusion as indicated in [o^, llOd l . The magnetic field in the upstream region can be 
large, and then the deflection angle of a particle after the shock crossing is large, too. Relativistic shock acceleration 
can operate in medium filled by pre- accelerated particles, and thus initial energy can be high. 

There are some other mechanisms of acceleration to e nerg ies up to ~ 10^^ eV relevant for AGN: unipolar induction 
and acceleration in strong electromagnetic waves (see [lOlj for description and references). The mechanisms of jet 
acceleration have very special status. 

The observed correlations between arrival directions of particles with energies (4 — 8) x 10^^ eV and BL Lacs [s^l 
imply the jet acceleration. This is because BL Lacs are AGN with jets directed towards us. For this correlation the 
propagation of particles (most probably protons) with energies above 4 x 10^^ eV must be (quasi)rectilinear, that can 
be realized in magnetic fields found in MHD simulations in [s^ (see however the simulations in [H, [13] with quite 
different results). 

An interesting mechanism of jet acceleration, called pinch acceleration, was suggested and developed in plasma 
physics [lO^. It is based on pinch instability well known in plasma physics, both theoretically and observationally. 
The electric current along jet produces the toroidal magnetic field which stabilizes the jet flow. The pinch instability is 
caused by squeezing the tube of flow by magnetic field of the current. It results in increasing the electric current density 
and magnetic field. The magnetic field compresses further the tube and thus instability develops. The acceleration 
of particles in [102| | is caused by hydrodynamical increase of velocity of the flow and by longitudinal electric fleld 
produced in the pinch. This process is illustrated by Fig.[THl The pinch acceleration has been developed for tokamaks 




FIG. 18: Pinch acceleration in the jet with electric current 



and was conflrmed there by observations. The generation spectrum is uniquely predicted as 

qgcn{E) cx with 7g = 1 + \/3 = 2.73. (21) 

Scaling a size of the laboratory tube to the cosmic jet, one obtains E'max exceeding 10^^ eV. The particle beam 
undergoes a few pinch occurrences during traveling along cosmic jet (a few kpc in case of AGN), and thus spectrum 
obtains a low energy cut at high value of i?min- This mechanism needs more careful study. 

We finalize this subsection concluding, that there are at least three possibilities for the broken generation spectrum 
with Qgan oc E~^ '^ dX E > Ec'- The acceleration by non-relativistic shocks with distribution of sources (more precisely 
emissivity) over i^max, the acceleration by relativistic shocks, where Ec = i?min naturally occurs due to the first 
u ^ d ^ u cycle as Tly^E^^ with E^^ being particle energy before acceleration, and in the pinch acceleration where 
7g = 1 + -s/S is rigorously predicted and i?min appears due to several pinch occurrences. The latter mechanism provides 
acceleration along a jet, necessary for correlations with BL Lacs. 



2. Spectra of UHECR from AGN 

We discuss here the diffuse energy spectra from AGN. They are model-dependent because one should specify the 
distance between sources d, the mode of propagation (e.g. rectilinear or diffusive) and the critical energy Ec in the 
broken spectrum of generation. The universal spectrum is valid when characteristic distance d between sources is 
smaller than other propagation lengths, most notably energy attenuation length Is^ttiE) and diffusion length Zdiff(-E). 
As was discussed in Section HVl the dip, seen at energies 1 x 10"'^* < < 4 x 10^^ eV, is the robust feature, which is 
not modified by any known phenomenon except presence of extragalactic nuclei with fraction higher than 20%. This 
fraction depends on model of acceleration and can be small for some acceleration mechanisms, e.g. for acceleration 
by relativistic shocks. 

The prediction of the GZK cutoff is model dependent: the shape of this feature can be modified by different values 
of Ejnax, by local overdensity /deficit of the sources (see Section fill Ap and by discreteness in the source distribution 
(see Fig. in]). It can also be mimicked by the shape of acceleration cutoff near E'max- 
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FIG. 19: Comparison of calculated AGN spectra with the data of AGASA and HiRes (left panel) and with data of Fly's Eye 
and Yakutsk (right panel). The generation spectrum is oc with Ec — 1 x 10^* eV and with m — 0. The propagation is 

diffusive with parameters as in Fig. [10] and with the Bohm diffusion in low-energy regime. 

The spectrum predicted for £■ < 1 x 10^® eV is also model dependent. At these energies for all distances d reasonable 
for AGN the universal spectrum is modified due to propagation in magnetic fields, because Idis is small. In Fig. [T7| 
(left panel) the spectrum for diffusion in random magnetic field is shown. At energies i? > 1 x 10^® eV the spectrum 
remains universal. 

We have to introduce in calculations the broken generation spectrum (|20p to provide the reasonable luminosities 
of AGN. Normahzed by the AGASA flux, the emissivity Co is given by 3.5 x lO'^^ ergs Mpc~'^yr^^, 1.6 x 10'''' ergs 
Mpc^^yr-i and 7.0 x 10'''^ ergs Mpc^^yr^^ for Ec equals to 1 x lO^^ eV, 1 x 10^^ eV and 1 x lO^*^ eV, respectively. 

To calculate the luminosities of AGN, Lp = Cq/us, we take the density of the sources from small-scale clustering 
as = 2 X 10"^ Mpc'^ [13, 111. The luminosities are found to be quite reasonable for AGN: 5.6 x 10'*^, 2.5 x lO**"* 
and 1.1 X 10**^ erg s~^ for Ec equals to 1 x 10^* eV, 1 x 10^'' eV and 1 x 10^^ eV, respectively. 

The calculated spectra are shown in Fig. [12] for propagation in magnetic field, as described in Fig. [TO] for the Bohm 
diffusion in the low-energy regime, for d = 50 Mpc, jg — 2.7, Ec — 1 x 10^^ eV and m = 0, inspired by BL Lacs, 
which show no positive evolution. The spectra for AGN evolutionary models are shown in Fig. [T3]in comparison with 
case m = 0. For calculations of spectra with different values of Ec see jl03j |. 



3. Transition to galactic cosmic rays 

For the diffusive propagation the transition is shown in Fig.[T7] In our calculations we followed Ref. [5^, considering 
the same case of diffusive propagation in random magnetic fields with basic coherent scale Ic — ^ Mpc and Bq = 1 nG 
on this scale. These parameters define the diffusion coefficient at very high energies, while at energies £^ < 1 x lO'^^ eV 
at interest we assume the Bohm diffusion. The distance between sources is fixed as d = 50 Mpc. As Fig. [17] shows 
the iron galactic spectrum given by curve "gal.Fe" describes well the KASCADE flux at i? « 1 x 10^^ eV and fits the 
Hall (drift) diffusion spectrum at £' > 1 x 10^^ eV. The predicted fraction of iron in the total flux is given in Table [T] 
as function of energy. This fraction is higher than in [5^ mostly because we used the HiRes data, instead of AGASA 
in 

The detailed study of transition for rectilinear propagation for different values of Ec has been performed in [l03j | . 
The fraction of iron- nuclei varies with Ec- 
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TABLE I: Fraction of iron nuclei in the total fiux as function of energy. 



E (eV) 


1 X 10^^ 


2 X 10" 


5 X 10" 


7 X 10" 


10i« 


Jfc/ Jtot 


0.97 


0.87 


0.49 


0.26 


0.04 



4- Signatures of AGN 

AGN are one of the best candidates for UHECR sources. They have high luminosities to provide the observed 
UHECR flux and have quite efficient mechanisms of acceleration. The predicted spectra are in good agreement with 
observations and transition to galactic cosmic rays is well described using diffusive or quasi-rectilinear propagation 
and broken generation spectrum. The density of sources ^ (1 — 3) x 10^^ Mpc"'^, found from small-scale clustering, 
corresponds well to space density of powerful AGN. However, in the diffuse spectrum we cannot predict any specific 
signature of AGN, and probably only observations of the direct flux of primary (e.g. protons) or secondary (e.g. 
neutrinos or gamma-rays) radiation can give the direct evidence for AGN as sources of UHECR. 

In this respect the correlations with BL Lacs [13, [13 is the direct indication to AGN as UHECR sources. The 
protons can be the signal carriers. The important point is that in structured (perforated) magnetic field, there are 
directions (holes) with weak magnetic fields in which BL Lacs can be seen [92] ■ In MHD simulation '35'| it was 
found indeed that in many large scale structures like voids and filaments the magnetic field is weak and protons can 
propagate there with small deflections. 

Correlation with BL Lacs imply the jet acceleration, because BL Lacs are AGN with jets directed to an observer. 
The pinch acceleration mechanism fits well the picture above, providing 7g = 2.7, high i?max and high £'minj needed 
for reasonable luminosity Lp. 

C. Gamma Ray Bursts 

GRBs are another candidates for UHECR sources. Like AGN they have high energy output and can accelerate 
particles to UHEs, as first have been proposed in [sol [ool. [9l| 

GRB starts with instant energy release Wtot, which occurs most probably due to SN explosion. This energy is 
injected into a small volume with initial radius Ri ~ 10^ cm in a form of a jet with small opening angle 6. The 
extremely large energy release within a small volume results in very high temperature of the gas, which consists 
mostly of photons and electron-positron pairs with a small admixture of initially injected baryons. The pressure of 
relativistic e"'"e~7 gas accelerates fireball, its temperature falls down and the Lorcntz factor rises. At later stages 
a fireball becomes the baryonic jet moving with Lorentz factor Fq = Wtot/MbC^ ^ 100, where Mi, is initial baryon 
loading. 

The jet is characterized according to observations jl04l | by small opening angle 9^4° and by corresponding solid 
angle fl = ttO'^. The external observer sees only part of the external fireball surface within the angle l/Fo < 9. Because 
if this fact, the convenient formalism is to consider the fictitious spherically-symmetric fireball with equivalent isotropic 
energy output Wtot and equivalent isotropic GRB rate vqrb (in units of Mpc"'^ yr~^)- The real energy output and 
GRB rate are Wtot = (f^/47r)Wtot and !>grb = (47r/r2)i^GRB, with the same emissivity C = WtotJ^GRB- 

According to observations the equivalent isotropic energy output Wgrb 1 x 10^"^ erg, and spherically-symmetric 
GRB rate is vgrb ~ 0.5 x 10~^ Mpc~^yr~^ with emissivity £grb ~ 0.5 x 10''^ erg Mpc^^yr^^. Th e more de tailed 
estimates differ considerably from this value: £gr,b « 0.6 x lO''^ erg Mpc^^yj.-! (j^^ Schmidt [lol [T06l} and 
^GR B ^ 1.3 X 10^'' erg Mpc~^yr~^ (Frail et al. |l04| '). Note, that difference in spectral intervals used in |l05l . [To^ 
and |l04l | cannot explain alone the differences between the above-cited emissivities. We conclude thus that £grb ~ 
(0.6 — 13) X 10^"^ erg Mpc'^^yr"^, to be compared with higher estimate by Waxman jlOTf £grb ~ (3 — 30) x 10"*^ erg 
Mpc~^yr~^. 

There are two sites of acceleration in a jet: mildly relativistic internal shocks inside a fireball [9ll| and ultra- 
relativistic external shock [90] . 

1. Internal shock acceleration and UHECR 
We shall describe here shortly the acceleration inside fireball, following [9l|, IIOTI . llOSj . 
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Gamma-radiation is produced as synchrotron emission of the shock-accelerated electrons at the stage when fireball 
moves with constant Lorenz-factor Fq = Wtot/MhC^. Acceleration occurs due to internal shocks produced by motion 
of fireball sub-shells with different Lorentz factors F. This phenomenon is assumed to be responsible for the observed 
GRBs and for their time substructure. In the observer frame the jet looks like a narrow shell with width A ~ R/^q, 
where i? is a radius, i.e. distance from outer surface to the source (the GRB engine). It is assumed that this shell 
consists of mildly relativistic sub-shells, moving relatively each other in the rest frame of the fireball. The inner sub- 
shells move faster than the outer ones, and their collisions produce mildly relativistic shocks. The GRB spikes arise due 
to shell substructure. The range of distances R-y where GRB is produced can be estimated from the observed minimum 
duration of spikes, Tgp > 1 ms, and duration of GRB, tgrb '--^ 1 — 10 s. Time between collisions of two sub-shells can 
be estimated in the observer frame [l08l | as 6t ^ S/v ^ TqS/c, where S is the distance between two colliding subshells. 
Thus, Rj = cSt is fimited between i?™'" ~ FgcT^p - 3 x 10" cm and i?^'^'^ - F^ctgrb 3 x lO^"* - 3 x lO^^ cm. The 
energetics is described by M^grb = /GReWtot, with /grb < 0.5 at least, because half of the total energy is released 
in afterglow. More realistically /grb ~ 0.2. 

The exit of all particles from the jet occurs only through the front spherical surface, provided b y cond ition 9 > I/Fq. 

Acceleration of electrons by internal shocks implies acceleration of protons. According to [9l|, Il07l | the generation 
spectrum cx E^^, -Bmax ~ 1 x 10^^ eV, E'min ~ TompC^ and Wp ^ W^grb (in the observer frame). 

As a matter of fact one must write the expressions above as £'min_2:' faXof^pf? a-nd Wp ~ /aW^GRB, where fa is a 
factor of adiabatic cooling. This phenomenon has been studied in 109*1 and adiabatic energy losses were found to be 
very severe. As was indicated there the efficient mechanism, which diminishes the influence of adiabatic cooling, is the 
neutron mechanism of UHECR exit from expan ding fireball. However, some restrictions obtained for this mechanism 
make its application rather limited. Waxman [107| fought back this criticism by argument that production of HE 
protons in the shell at later periods corresponding to a distance R^^^ eliminates the problem of adiabatic cooling. 
We shall argue below that this problem exists. 

There are three stages of adiabatic energy losses. The first one occurs during acceleration process and it operates at 
distances R™™ ^ R ^ i?™'*'^. The adiabatic energy losses are described in the standard shock-acceleration equation by 
the term (p/S)div u df{r, t,p)/dp, where u is hydrodynamical velocity of the gas flow and /(r, t^p) is the distribution 
function. In our case the adiabatic energy losses are large due to expansion of the shell. Adiabatic energy losses 
diminish the fraction of the shock energy transferred to accelerated particles, which results in relation Wp = /i^^W^GRB, 

where /i^' < 1 is adiabatic factor. Qualitatively it can be explained in the following way. The first portion of 
accelerated particles would suffer severe adiabatic energy losses fa ~ R^^'^/R^'^^, if further acceleration were absent. 
In fact, these particles are re-accelerating and loosing energy simultaneously, as corresponding diffusion equation 
describes it. This process diminishes only the fraction of energy fa transferred to accelerated particles. It can be 
estimated as fa ~ R'^^/R^'^^, but for accurate estimate the solution of diffusion equation is needed. Unfortunately, 
this model exists only in the form of estimates. 

The second stage starts after acceleration ceases. The end of GRB radiation signals this moment. Magnetic field 
in the expanding shell slowly diminishes, accelerated electrons radiate away their energy fast, but protons are still 
confined in the shell, forming relativistic gas and loosing energy adiabatically. Duration of this stage is different for 
protons of different energies and the spectrum oc E~'^ is distorted. The protons with lO'^^ — 10^^ eV are the first 
to accomplish this stage and enter the regime of free expansion (stage 3). This energy range is most interesting for 

(2) 

UHECR . The adiabatic factor fa depends on the time of diminishing of equipartition magnetic field. The detailed 
calculations are needed to evaluate /i^'' for the range 10^^ — 10^^ eV. 

At the third stage when the shell of the protons with energies 10^^ — 10^^ eV propagates quasi-rectilinearly in the 
external magnetic field, the energy losses are also present. They are caused by the collective effect of transferring 
proton momenta to the macroscopic volume of gas due to scattering in magnetic field. 

As a toy model, let us consider the spherically symmetric shell of protons with energies 10^^ — 10^^ eV, and make an 
order of magnitude estimates. At the distance of order of gyroradius r^, a proton transfers half of its initial (radial) 
momentum to the gas. Since protons propagate in the radial direction almost with speed of light, it follows from 
causality that volume of gas, which absorbs the lost momentum cannot be larger than ~ r^. The equation of radial 
momentum balance reads 

0.5Wp/c ~ Trlpgc, (22) 

where pg ~ 10^^^ g/cm'^ is the density of the ambient gas and F is Lorentz factor, obtained by this gas volume. For 
the spectrum c>c E^^ the total energy of protons with 10^^ < E < 10^^ eV is Wp « 0.2iyGRB- For magnetic field in 
th e pre supernova wind at dis tanc e r ~ 10^^ cm we accept the value of magnetic field i?o ~ 1 G, following the estimate 



of [100| 0.1 - 10 G (see also [l07|). Then for W^grb ^ 1 x lO^^ gj.g we obtain F - 10^ and for Wqrb ~ 5 x lO^o erg 
(jet value) F ~ 10^. In fact for the jet case the effective mass of the target is much less and F is larger. Such gas shell 
produces the relativistic shock and considered scenario radically changes. 



26 



If considered UHE protons capture electrons from the media with its frozen magnetic field, the pressure of the beam 
to the gas further increases. 

In fact, we considered above only one mechanism for energy transfer to the gas. It is clear, that nature of this 
phenomenon is very general: the energy density of proton beam is so large, that part of its energy is transferred to 
magnetized ambient plasma. Only when this density drops to some critical value and collective effects disappear, the 
protons can propagate as individual particles, not exciting the ambient plasma. 

Therefore, the total adiabatic factor fa cannot be of order of 1. 

In the left panel of Fig. [^D] the calculated flux of UHE protons accelerated by internal shocks is presented by curve 1. 
It corresponds to generation spectrum cx with i?max = 1 x 10^^ eV and i?min = FompC^ = 1 x 10^^ eV. This is the 
case of absence of adiabatic energy losses /a = 1. The spectrum is calculated in the standard way taking into account 
the energy losses on CMB radiation. The emissivity needed to fit the observational data is Ccr = 2.2 x lO'^^ erg 
Mpc~'^yr~^. It must be equal to Wpi^cRB and thus it should not exceed £grb given above. As a matter of fact it is 
by factor 20 - 300 higher. If adiabatic energy losses for protons with energies 10^^ — 10^"'^ eV is /a(> 1 x 10^^), the 
predicted fiux in this energy range is additionally suppressed by this factor and discrepancy increases correspondingly. 

If to assume fa{E < 1 x 10"'^^ eV) <^ 1, while fa{E > 1 x 10^^ eV) « 1, the energy balance is not changed, because 
the total energy (including Wp and that adiabatically lost) remains the same. 




E, eV E, eV 

FIG. 20: UHECR from GRBs. In the left panel are shown the fluxes produced by internal shocks (curve 1) with generation 
spectrum oc E~^, JSmax ~ 1 x lO^'^ eV and -Emin ~ 1 x 10^^ eV and by external relativistic shocks (curve 2) with generation 
spectrum oc -Emax ~ 1 x 10^^ eV and -Emin ~ 1 x 10^'' eV. In the right panel the spectra are shown with the assumed 

distribution of sources over Emax, like in the case of AGN. The effective generation spectrum ai E > Ec is assumed to be 
oc E-^-'' with Ecr^lx 10^* eV. 

In the right panel of Fig. [201 we show the spectrum with the assumed distribution of sources over i?max like in case 
of AGN. In this case the required emissivity is higher Cqr — 3.5 x lOf® erg Mpc~'^yr~'^. 

These results have been already presented in [7l|. Waxman [107^ did not rmderstand that we considered two 
different cases (now being put separately in two different panels in Fig. I20p and argued in fact only against the case 
shown in the right panel of Fig. 1201 The left panel shows the same case as considered by Waxman. Our discrepancy in 
emissivity is explained as follows. We normalize our calculations by emissivity Cp (the energy released in cosmic ray 
protons per unit comoving volume and per unit time) a nd c ompare it with >Cgrb — W^GRBi^GRB- Waxman normalizes 
his flux by value E^hcR/dE — 0.8 x lO^'' erg/Mpc^yr jl07l |. To obtain emissivity, this value should be multiplied by 
ln(£^„iax/-Emin), which is equal to 23 for E^max = 1 x 10^-^ eV and i?min = 1 x 10^^ eV. The obtained UH ECR emissivity 
roughly agrees with our value, but it exceeds by an order of magnitude what Waxman cites in jl07l | as the highest 
value of GRB emissivity (3 x 10'*'* erg Mpc~'^yr~*). Note that the comparison is made for unrealistic case of absence 
of adiabatic energy losses. 
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2. External shock acceleration and UHECR 

The flow of ultra-relativistic baryonic jet is terminated by the external ultra-relativistic shock propagating in pre- 
supernova wind. The acceleration by such shocks has been described (with relevant references) in subsection I VII BT] 
As a shock propagates in interstellar medium, the protons, experiencing the different number oi u ^ d ^ u cycles, 
move together with the shock ahead of it (upstream) or behind it (downstream). The "fresh" particles, which undergo 
one or a few (n) cycles have energy E ^ T^^2"-Ei, the "old" particles approach ii^max ~ sBRgiiTsh [SJ]' where Ei is 
initial proton energy upstream, i?sh is a radius of the shock (the distance from GRB origin), and Fgh 100 is the 
Lorentz-factor of the shock. The acceleration time ta ^ E / (ecSFgh) .94] is smaher than the age of the shock t ~ Rsh/c 
for all particles except those with E ~ i?max- 

At this stage of our study we are interested only in the total energy transferred to accelerated protons, disregarding 
what happens to these particles later. In other words we are interested only in energy balance, using the common 
assumption Wp ~ IVgrb for the total energy Wp transferred to accelerated protons. 

For the typical afterglow distance Rsh ~ 10^® cm, when the shock moves with constant Lorentz factor Fgh ~ 100, we 
shall use the following parameters for the upstream HE protons: the power-law spectrum c>c E^'^n with the "standard" 
(see subsection IVII B jg = 2.2, with conservative £^min ^ ^sh^i ~ ^^^^ '^i*'^ -^'max = eBRs^Tsi^/^/S w 

1 X 10^-'^ eV, valid for very large i? w 5 G. The fraction of the total energy Wp in the form of particles with 
Ep > 1 X 10^^ eV is /uHECR — 0.039. In exotic case of existence of pre-accelerated particles in pre-supernova wind 
with Ei ^ 10^^ eV, the minimal energy is i?min ~ 10^^ eV and /uhecr = 0.29. 

As a matter of fact, part of accelerated particles are located downstream. Because they are isotropized in the rest 
frame, the minimum energy in the observer frame is shifted to 

i?mi„ = ^^(l + F;VF,\), (23) 

where F^ is the Lorentz factor of proton in the rest system. The energy (|23p corresponds to a proton moving backward 
relative to shock. Additionally the number of particles with maximum energy diminishes due to angular distribution. 
It results in diminishing of /uhecr! we shall not take this effect into account. 

There is one more effect which diminishes /uhecr for downstream protons. The highest energy protons can scatter 
to backward direction and escape through the rear surface of the shell, while the low-energy protons are better confined 
in the downstream region. In the laboratory frame the backward moving protons have energy about 1 GeV, as given 
by Eq. The difference of proton energy (in the laboratory frame) is transferred to the shell in the process of 

magnetic scattering. This energy will go back to particle acceleration, and most of it will be taken by low-energy 
particles. 

Adiabatic energy losses, not taken into account in the relativistic shock acceleration, are discussed in previous 
subsection. 

In realistic models the shock propagates in diminishing magnetic field, e.g. B{R) = BqRq/R [lOOi] . valid for the 
stellar wind. For such field it is easy to calculate the deflection angle 9 for proton propagating from i?sh to R, 

e«^%^lni?/i?sh, (24) 
III 

and the deflection angle at which the shock overtakes the proton, 9c — VS/Fsh, is the same as for constant magnetic 
field. The maximum acceleration energy i?max determined from the condition 9 = 9c remains also approximately the 
same. The exit of particles in the interstellar medium occurs when Tgh diminishes, most notably at the Sedov phase. 

In Fig. [201 the calculated spectrum of UHECR from acceleration by external shock is given by curve 2 for 7^ = 2.2 
and E'max = 1 x 10^^ eV. The UHECR emissivity needed to fit the flux at E' = 1 x 10^^ eV to the observational data 
is C{> 1 X lO^^eV) = 4.7 x 10"*^ erg Mpc"'^ yr~^. The lower limit to HE proton emissivity (without adiabatic energy 
losses and other effects discussed above) is given by 

Cp > £(> 1 X 10i9eV)//uHECR = 1.2 x 10*^ ergMpc'^r"' (25) 

for £^„iin = 1 X 10^^ eV, and Cp > 1.6 x lO''^ erg Mpc~^ yr~^ for exotic case with i?min = 1 x 10^'' eV. The emissivity 
(1^5)1 is 2 - 3 orders of magnitudes higher than the observed GRB emissivity (0.6 — 13) x 10'*^ erg Mpc^^yr^^. The 
exotic assumption of -Bmin = 1 x 10^^ eV reduces this discrepancy only by one order of magnitude. 

3. Discussion and conclusion 

The calculations above agree with earlier works [7l|, [t^, IllOl llll| (note the emissivities given in the last work of 
this list). 
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The main diff eren ce between calculations [7l|, [tI, IllOl on the one hand, and those by Waxman (e.g. |l07| ') 

and Vietri (e.g. [lOClj |l. on the other hand, consists in the following. 

While Waxman and Vietri compare GRB emissivity with emissivity of the observed UHECR (Vietri a,t E > 
1 X 10^^ eV), we compare it with emissivity, which corresponds to energy transferred to all protons in the process of 
their acceleration. Apart from spectral fraction of energy, it includes also the energy lost by protons before they escape 
to extragalactic space (adiabatic energy losses, energy transferred to the magnetized plasma outside the jet and other s 
discussed above). Waxman for UHECR energy input uses the value E-^dficR/dE = 0.8 x 10^'' erg Mpc~'^yr~^ jlOTj . 
which after integration from E'min = 1 x 10^^ eV to -Emax = 1 x 10^^ eV corresponds to emissivity Cp = 1.8 x 10'*^ erg 
Mpc~^yr~^ to be compared with the observed GRB emissivity (0.6— 13) x lO''^ erg Mpc~^yr~^. Because of adiabatic 
energy losses, £„ is expected to be in fact much higher than the above-mentioned value. 

Vietri et al. |lOOl | indicate UHECR energy release required by observations as 1 x 10^^ erg Mpc^'^yr^^ for energy 
interval 1 x 10"'^'^— 1 x 10^^ eV (compared with our value 1.6 x lO''^ erg Mpc~^yr~^ for the same energy interval). If one 
includes the energy input to all accelerated protons, which accompany ultra-relativistic shock, from Emin ~ ^sh-^i 
Emax, this emissivity increases at least by factor 7.3, in conflict with observations by factor 60. There is some difference 
with our work in numerical calculations: Vietri et al. uses the HiRes data, we - energy-shifted data by AG AS A 
and HiRes, which are in good intrinsic agreement (see Fig. 120^ : Vietri et al. introduce corrections due interaction 
fluctuations, we do not need to do it, because we normalize our calculations by observational data at £' = 2 x 10^^ eV, 
where fluctuations are absent. The maximum acceleration energy i?max = 1 x 10^"'^ eV at Rsh = 1 x 10^^ eV requires 
according to our calculations (it is by factor \/3 less than in [94]) magnetic field in the p resuperno va wind at this 
distance ~ 6 G, which is very large and marginally possible only for the SupraNova model jll2llll3j . 

Our calculations require the UHECR emissivity £(> 1 x lO^^eV) = 4.4 x 10^^ erg Mpc~'^yr~^ for internal shock 
acceleration (7g = 2.0 and £^max = 1 x 10^^ eV) and 4.7 x 10^** erg Mpc~^yr^^ for external ultra-relativistic shock 
acceleration (7^ = 2.2 and i^max = 1 x 10^^ eV). Including the energy injected in low energy protons (not necessarily 
coming outside), this emissivity should be a minimum of 10 - 30 times higher, and thus contradiction with GRB 
emissivity reaches conservatively the factor 30 - 100. A natural solution to this problem is given by an assumption 
that accelerated protons carry ~ 100 times more energy than electrons, whose radiation is responsible for GRBs. 
However, for 4° beaming it increases the required energy output of SN as GRB source from VFsn ~ 5 x 10^" erg to 
l()52 _ ^q53 gj-g^ which creates the serious problem for SN models of GRBs. This problem exists even for the case 
when one takes into account only "observed" UHECR emissivity £(> 1 x lO^^eV). Then the minimal energy balance 
for internal shock model includes £grb = 1-3 x lO**^ erg Mpc~^yr~^, /^afterglow ~ 'Cgrb, 1 x lO^^eV), and 

Ci, w £(> 1 X lO^^eV), which results in SN energy output 4 x 10^^ erg. For external shock model the energy balance 
does not include neutrinos and it results in Wsn ~ 3 x 10^^ erg. 

VIII. CONCLUSIONS 

This work is naturally subdivided into two major parts: (i) the model- independent analysis of spectral signatures of 
UHE proton interaction with CMB, and (ii) the model-dependent analysis of transition from extragalactic to galactic 
cosmic rays and model-dependent study of most probable UHECR sources: AGN and GRBs. 

In the first part the number of assumptions is minimal: we assume the power-law generation spectrum of extra- 
galactic protons oc E~'^^ and analyze uniform distribution of the sources with a possible distortion in the form of 
large-scale source inhomogeneities and local overdensity or deficit of the sources. The reference spectrum is given by 
the universal spectrum, calculated for the homogeneous distribution of the sources, when separation of sources c? — > 0. 
According to the propagation theorem [5^, in this case the spectrum does not depend on propagation mode, e.g. it 
is the same for rectilinear and diffusive propagation. We analyze also the spectra with different source separations, 
giving emphasis to separations d ~ 30 — 50 Mpc, favorable by clustering of UHECR [2^ [27l. l28t . 

Calculation of the universal spectrum is based on equation of conservation of particles with continuous energy losses 
of protons interacting with CMB. We performed the new calculations of energy losses (first presented in [lH), which 
agree well with (4ll \5T\ . see Fig. [T] The technical element needed for calculations, the ratio of energy intervals at 
generation and observation dEg/dE, is calculated by the new methods (see Appendix iBjl . The analysis of spectral 
features is convenient to perform with help of modification factor, r](E) = Jp(E)/ Jp'^"^{E), where the spectrum Jp{E) 
is calculated with all energy losses included, and unmodified spectrum Jp'""(i?) - only with account of adiabatic 
energy losses (red shift). 

There are four major signatures of UHE protons interacting with CMB: GZK cutoff, bump, dip and the second dip. 

GZK cutoff [8, 9] is the most prominent signature, which consists in the steepening of the proton diffuse spectrum 
due to photopion production. The beginning of the GZK cutoff at i?GZK ~ 4 x 10^^ eV is difficult to observe. The 
quantitative characteristic of the GZK cutoff is given by energy E1/2, where the flux with the cutoff becomes lower 
by factor 2 in comparison with power-law extrapolation from lower energies. This quantity is possible to observe 
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only in integral spectrum. Practically independent of 7g i?i/2 = 5.3 x 10 eV. Data of Yakutsk and HiRes are 
compatible with this value (see Fig. 2]), while data of AG AS A contradict it. As to spectral shape of the GZK cutoff, 
we found it very model-dependent: it depends on separation of sources d, on local source overdensity or deficit, on 
maximum acceleration energy -Emax, and being so uncertain, the GZK steepening can be imitated by -Emax acceleration 
steepening. 

Bump is produced by pile-up protons, experiencing the GZK process. These protons do not disappear, they only 
shift in energy towards that where the GZK cutoff begins. The bump is very clearly seen in the spectra of individual 
sources (see left panel of Fig. [7]), but they disappear in the diffuse spectrum (right panel of Fig. [7]), because individual 
bumps are located at different energies. 

Dip is the most robust signature of UHE protons interacting with CMB. It has very specific shape (see Fig. [2]), 
which is difficult to imitate by other processes, unless one has a model with many free parameters. Dip is produced by 
pair production on the CMB photons p + j—^p + e^ + e^ and is located in energy interval 1 x 10^^ — 4 x 10^^ eV. Dip 
has two flattenings (see Fig. [2]), one at energy E ^ 1 x 10^^ eV, which automatically reproduces the ankle, well seen 
in the observational data (see Fig. [8]), and the other at i? ^ 1 x 10^^ eV, which provides transition from extragalactic 
to galactic cosmic rays. The prediction of the dip shape is robust. It is not noticeably modified by many phenomena 
included in calculations: by different distances between sources (Fig. [5]), by changing the rectilinear propagation of 
protons to the diffusive (Fig. [TO)) , by different i?max (Fig- E]) , by local overdensity and deficit of the sources (Fig. [6]) , 
by large-scale source inhomogeneities and fluctuations in p7- interactions fFig. [T5)) . 

The dip is very well confirmed by data of Akeno-AGASA, HiRes, Yakutsk and Fly's Eye (see Fig. [inifor the latter), 
while data of Auger at this stage do not contradict the dip (Fig. [5]). The agreement with the data of each experiment 
Akeno-AGASA, HiRes and Yakutsk is characterized by « 19 — 20 for about 20 energy bins, and with only two free 
parameters, jg and the flux normalization constant. The best fit is reached at 7g — 2.7, with the allowed range 2.55 
- 2.75. 

The dip is used for energy calibration of the detectors, whose systematic energy errors are up to 20%. For each 
detector independently the energy is shifted by factor A to reach the minimum x^- We found AAg = 0.9, Xm — 1-2 and 
Ava = 0.75 for AGASA, HiRes and Yakutsk detectors, respectively. Remarkably, that after this energy shift the fluxes 
and spectra of all three detectors agree perfectly, with discrepancy between AGASA and HiRes at i? > 1 x lO^'' cV 
being not statistically significant (see Fig. fTTj) . The AGASA excess over predicted GZK cutoff might have the statistical 
origin combined with systematic errors in energy determination, and the GZK cutoff may exist. 

The difference in energy shifts between fluorescent energy measurement (HiRes) and on-ground measurements 
(AGASA, Yakutsk) with A > 1 for HiRes and A < 1 for AGASA and Yakutsk may signal a systematic difference in 
energy determination by these two methods. 

The excellent agreement of the dip with observations should be considered as confirmation of UHE proton interaction 
with CMB. 

For astrophysical sources of UHECR the presence of nuclei in primary radiation is unavoidable. Whatever a 
source is, the gas in it must contain at least hehum of cosmological origin with ratio of densities riHe/'^H = 0.079. Any 
acceleration mechanism operating in this gas must accelerate helium nuclei, if they are ionized. The fraction of helium 
in the primary radiation depends critically on the mechanism of injection. Fig. 1131 shows that while 10% mixing of 
helium with hydrogen in the primary radiation leaves good agreement of dip with observations, 20% mixing upsets 
this agreement. One may conclude that good agreement of dip with observations indicates a particular acceleration 
mechanism, more specifically - an injection mechanism. One may also note that the shape of the dip contains the 
information about chemical composition of primary radiation (see Fig. [T^ . 

The low-energy flattening of the dip is seen for both rectilinear and diffusive propagation (see Fig. [T6)) . It is 
explained by transition from adiabatic energy losses to that due to e+e^-pair production. Low-energy flattening of 
the extragalactic spectrum provides transition to more steep galactic component, as one can see it in the left panel 
of Fig. [T71 Note, that when spectrum is multiplied to E^-^ the flat spectrum looks like one raising with energy. In 
Fig. [17] the transition is shown for the diffusive spectrum. This model of transition implies that at energy below 
Eti- « 5 X 10^^ eV the galactic iron is the dominant component, while at the higher energy extragalactic protons with 
some admixture of extragalactic helium dominate. The observational spectrum feature, which corresponds to this 
transition, is the second knee. The prediction of this model, the dominance of proton component at i? > 1 x 10^^ eV, 
is confirmed by data of HiRes [l^, [3, HiRes-MIA [l^, by Yakutsk [l^ and does not contradict to the Haverah 
Park data at > (1 — 2) x 10^^ eV [if. However, the AkenoflTt and Fly's Eye Q data favor the mixed composition 
dominated by heavy nuclei. 

The second dip is produced by interplay between pair production and photopion production. It is a narrow feature 
at energy £'2dip = 6.3 x 10^^ eV, which can be observed by detectors with good energy resolution. If detected, it gives 
the precise mark for energy calibration of a detector. 

The UHECR sources have to satisfy two conditions: they must be very powerful and must accelerate particles to 
Emeix > 1 X 10^^ eV. There is one more restriction [H, [H, [23, HE HI, [13, El] , coming from observation of small-scale 
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clustering: the space density of the sources should be (1 — 3) x 10~^ Mpc~^, probably with noticeable uncertainty 
in this value. Thus, these sources are more rare, than typical representatives of AGN, the Seyfert galaxies, whose 
space density is ~ 3 x 10^'' Mpc^'^. The sources could be the rare types of AGN, and indeed the analysis of (30l | 
show statistically significant correlation between directions of particles with energies (4 — 8) x 10"'^^ eV and directions 
to AGN of the special type - BL Lacs. The acceleration in AGN can provide the maximum energy of acceleration 
to i?max ~ 10^^ eV. The appropriate mechanisms can be acceleration by relativistic and ultra-relativistic shocks (see 
section rVlIB 1[) . These mechanisms typically provide spectra flatter than one with 7g « 2.6 — 2.8, needed for agreement 
of the dip with data. The distribution of the sources over Emax suggested in [2^ solve this problem: starting from 
some energy Ec the spectrum becomes steeper. In this case the diffuse spectrum below Ec has canonical generation 
index jg = 2.0 or = 2.2, and above this energy 7^ « 2.6 — 2.8 to fit the observations. The energy Ec is a free 
parameter of the models, and we keep it typically as 10^^ — 10^* eV. Such broken generation spectrum relaxes the 
requirements for source luminosities. For example, for 7^ = 2.7 and Ec = Ix 10^® eV, the typical emissivity needed to 
fit the observed flux is £0 = ngLp ~ 3 x 10**^ erg Mpc~^yr~^, which is held for the source density = 2 x 10~^ Mpc~'^, 
estimated from the small-angle clustering [l^, [l^, and for the source luminosity Lp ~ 6 x lO'^^ erg/s, very reasonable 
for AGN. The calculated spectra are in a good agreement with observations (see Fig.[T9|). The transition from galactic 
to extragalactic cosmic rays in this model occur at the second knee (see left panel of Fig. [T7)) . 

Another potential UHECR sources are GRBs. They hav e lar ge energy output and can accelerate particl es t o 
UHEs. Acceleration occurs at the multiple inner shocks (oil. Ill4| and at ultra-relativistic external shock [ool Il00j |. 
The maximum acceleration energy can be as high as ~ 1 x 10^^ eV. In case of ultra-relativistic shock with Lorentz 
factor of shock Tgh ~ 100 and at distance Rsh ~ 1 x 10^^ cm, typical for afterglow region, the maximum energy 
£^max = (l/V^)eBi?shrsh reaches 1 x lO^i eV in case of very strong external magnetic field i? ~ 6 G. Such strong 
magnetic field can be margi nally pro vided by two combined effects: strong field in presupernova wind (possible only 
in the model of SupraNova |ll2l . IllSj ) combined with preshock amplification of the field. 

T he en ergy output of GRBs necessary to provide the observed UHEGR flux imposes the serious problem [7l|, I?!, 
, llll| . It is usually assumed that energy output in the form of accelerated protons is of order of magnitude of that 
observed in GRB photons. The observed emissivity of GRBs is 

C - (0.6 - 13) X 10*3 ergMpc-Vr"\ (26) 

where the lower value is from works by Schmidt |l05l . Il06j | and the large value - from work by Frail et al. |l04{ . Both 
authors used in their estimates the similar energy interval of photons, (10 - 1000 keV) and (0.2 - 2000 keV). 

Our calculations of UHECR spectra for internal shocks (see Fig curve 1) with 7^ = 2.0 and E'max = 1 x 10^^ eV 
give UHECR emissivity £(> 1 x 10^^ eV) = 4.4 x 10'^'' ergMpc"^yr-\ by factor 3.4 - 73 higher than If to 

include ~ -Cuhecr contradiction with (|26p increases to 6.8 - 150. If to include in balance the protons with energies 
down to E'min '^-^ 1 X 10^^ eV (not necessarily coming out) the contradiction increases by another factor 5. Adiabatic 
energy losses further increase this discrepancy. 

Our calculations for external shock acceleration (see curve 2 in Fig 1^0)) with 7g = 2.2 and i?max = 1 x 10^^ eV 
results in UHECR emissivity £(> 1 x 10^^ eV) = 4.7 x lO^'^ ergMpc"^yr-\ i.e. by factor 3.6 - 78 times higher than 
(p6)) . It increases by factor 26, if to include protons with energies down to i^min ^ ^sh-^i ^^^^ by factor 3.5 

for unrealistic case of E'min ~ 1 x 10^^ eV. Adiabatic energy losses are less than in case of internal shock acceleration. 

This energy crisis can be resolved, assuming that accelerated protons in GRBs take away ^ 100 times more energy 
than electrons responsible for GRB photons. However, this assumption dramatically affects the status of SN origin 
for GRBs , even in the case of narrow 4° beaming. 

A. Predictions for the Auger measurements 

What do we predict for the Auger observations on the basis of our analysis? 

• We predict the dip at energy 1 x 10^® — 4 x 10^^ eV not only as the spectrum shape, but with the absolute 
values of flux, based on the agreement of AGASA-HiRes- Yakutsk fluxes after energy calibration by the dip (see 
Fig. [TT]) . This spectrum is displayed in the left panel of Fig. [5TJ The distortion of the shape of the predicted 
spectrum can occur only due to presence of extragalactic nuclei. 

• The above-mentioned distortion is important as a tool for measuring the chemical composition of UHECR. The 
chemical composition is the most difficult problem for UHECR experiments, in which little progress has been 
reached during last 30 years. The measurement of the dip shape gives the precise method of determination of 
10 - 15 % of admixture of helium in primary UHECR (see Fig. [T2)) . The distortion of the proton dip most 
noticeably occurs at energy i? > 3 x 10^^ eV, (see Fig. [T3]) where precise measurements of Auger are expected. 
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FIG. 21: Predictions for the Auger detector. In the left panel the diffuse spectrum prediction is displayed for different distances 
d between sources, which provide the largest uncertainties. The spectrum is normalized by the energy-shifted AGASA data 
which coincides well with HiRes and Yakutsk data (see Fig. Ilip . In energy range 1 x 10^* — 7 x 10^^ eV the uncertainties in 
the predictions are small and Auger should observe the beginning of the GZK cutoff at i5 > 5 x 10^® eV. In the right panel the 
spectra with bumps from the nearby sources are displayed (7^ is the generation index). 

• Normalization of energies measured by AGASA, Yakutsk and HiRes detectors with the help of the dip (see 
section IIV CI) implies that energy measured by fluorescent method is higher than that by on-ground method 
(the energy shift for AGASA and Yakutsk is A < 1, while for HiRes A > 1). With help of the hybrid events 
Auger can resolve this problem. 

• The second dip (see Fig. I24p may be marginally observed and used for energy calibration of the detector. 

• Yakutsk and HiRes data give evidence for the numerical characteristic of the GZK cutoff E1/2 — 5.3 x 10^^ eV 
in the integral spectrum (see Fig. 0]). In the Auger data this value can be measured more reliably. 

• The predicted shape of the GZK steepening has many uncertainties (see section UlI Al for the discussion). How- 
ever, beginning of the GZK cutoff in the narrow energy interval (4—7) x 10^^ eV is predicted with the high 
accuracy (see left panel in Fig. l2ip . The precise spectrum measurement in this energy range can give the decisive 
proof of the GZK cutoff. 

• The mass composition at 1 x 10^^ < £' < 1 x 10^^ eV measured by the fluorescent technique can distinguish 
between models of the second knee transition from galactic to extragalactic cosmic rays (protons with 10 - 20% 
of helium) and the ankle transition (iron dominance) 

• The bumps in the spectrum of single sources, if detected, can be be observed (see Fig. [5T|). 
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APPENDIX A: CALCULATIONS OF ENERGY LOSSES 

We consider here some details of calculations of energy-losses due to production of e+e~-pairs and pions. 
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Pair production energy loss of ultrahigh-energy protons in low-temperature photon gas, e.g. CMB, 



P + ICMB -*p + e^ + e~ 



(Al) 



has been previously discussed in many papers. The differential cros s-sec tion for this proc ess in the first Born ap- 
proximatio n wa s originally calculated in 1934 by Bethe and Heitler |ll5l | and Racah [ll6j |. In 1948 Feenberg and 
Primakoff |ll7l | obtained the pair production energy loss rate using the extreme relativistic approximation for the 
differential cross-section. In 1970 the accurate calculation was performed by Blumenthal [1^. Later some analytical 
approximations to differential cross-sections were obtained in Ref . jll8l | . 

All authors neglected the recoil energy of proton, putting rup — > oo, the effect being suppressed by a factor of 
nie/mp « 5 X 10~^. 

In spite of the fact that all calculations actually used the same Blumenthal approach, there are noticeable discrep- 
ancies in the results of different authors; they are displayed in Fig. lb of the Ref. [4^ . 

To clarify the situation we recalculate the pair production energy loss of high-energy proton in the low-temperature 
photon gas. In contrast to Ref. [35] we use the first Born approximation approach of Ref. ^40.] taking into account the 
finite proton mass. The exact non-relativistic threshold formula with corrections due to different Coulomb interactions 
of electron and positron with the proton (see e.g. Ref. (119]) is used. No series expansions of cr(er), where Er is the 
photon energy in the proton rest system, are involved in our calculations. The average fraction of proton energy loss 
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FIG. 22: In the left panel is shown the average fraction of the incident proton energy Ep carried away by e e~ pair as a 
function of the photon energy Sr in the proton rest frame. The right panel displays the product of fraction of energy lost, 
(1 — x), and the cross-section for pair production, 7p — » e'^e~p, or photopion production p -\- ^ —^ X as function of Er- 



in one collision with a photon is defined as 

/ dx(l-x)^%^, (A2) 
(t(£,.) J ax 

where x — E'p/Ep, with Ep and Ep being the incident and final proton energies, respectively, in the laboratory system. 
The upper and lower limits on fractions x are given by 

a;max(e^) = Xcisr) ± ydsr), (A3) 



1 / m'i — ml \ / 

Xc{Sr)^^[l+ ; yc{Ee)^Jxl{Sr)-^-f • (A4) 

2 \ + zmpEr J y m,^ + 2mper 

Our strategy was to calculate {1 — x){er) by performing the direct fourfold integration of the exact matrix element 
over the phase space. It should be noted, that direct numerical integration, especially at high energies, is difficult in 
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this case because of forward-backward peaks in the electron-positron angular distributions. To overcome this problem, 
we performed two integrations over polar and azimuth angles in the e^e~ subsystem analytically. This was facilitated 
by using of the MATHEMATICA 4 code. The residual two integrations over energy and scattering angle in the initial 
P7 subsystem were carried out numerically. We calculate simultaneously the total cross-section for pair production. 
The accuracy of our calculations was thus controlled by comparison of the calculated total cross-section with the 
well-known Bethe-Heitler cross-section. 

The average fraction of proton energy lost in one collision with a photon is plotted in Fig. [51] as a function of the 
photon energy in the proton rest system. 

The product of this fraction and the total cross-section for pair production is shown in Fig. [221 This function should 
be integrated over the photon spectrum to obtain the average energy loss due to pair production in the photon gas 
with this spectrum. 

The comparison of our calculations with Ref. 42] shows the negligible difference (see Fig.[T]). 

Calculations of photopion energy losses are described in section Hi Al In the right panel of Fig. [22| we show {1 — x)a 
for these energy losses as function of Sr- 



APPENDIX B: CONNECTION BETWEEN ENERGY INTERVALS AT EPOCHS OF PRODUCTION 

AND OBSERVATION 

If we consider the protons with energy E in the interval dE at the epoch with red-shift z — 0, what will be the 
corresponding interval of generation dEg at epoch z, when energy of a proton was Eg{z)7 The connection between 
these two intervals is given by Eq. (36) of Ref. [4lj. Here we shall confirm this formula using a different, more simple 
derivation. Note, that intermediate formulae (40) and (41) used for derivation of final Eq. (36) in Ref. [4l| have a 
misprint: the correct power of (1 -I- z) term there is 3, not 2. 

Regarding the proton energy losses of a proton on CMB at arbitrary epoch with red-shift z, we shall use, as in 
Section [ira the notation b{E,z) = -dE/dt and P{E,z) = -{l/E)dE/dt with b{E,z) ^ {I + z)^6o[(l + z)E]. Here 
and henceforth ho{E) and (3q{E) are used for energy losses at z = 0. The energy loss due to red-shift at the epoch z 
and the Hubble parameter H{z) are given by 



{dEldt\_^^ = -EH{z), H{z) = Ho^n^il + zY + n^. (Bl) 
The energy of a proton at epoch z, 

Eg{z)^E + j\t [{dEldt\_,^ + {dE/dt)^^^] , (B2) 

can be easily rearranged as 

Eg{z)^E+ rJ^Egiz')+ r dz'\±4rbo[{l + z')Egiz')], 
Jo ^ + z Jo J^[z ) 

using dt/dz = -H-^{z)/{l + z) and Eq. (O. 

Differentiating Eq. (|B3p with respect to E, one finds for energy interval dilation y(z) = dEg{z) / dE: 

y(^) = 1 + X YT^yi' ) + ldz -j^y{z ) [-^ ] ■ (B4) 



(B3) 



E' = (l+z')Eg(z') 



Corresponding differential equation is 

1 dy{z) _ 1 (1 + z)2 fdbo{E') 



y{z) dz l + z H{z) V dE' 

The solution of Eq. (|B5|) is 



E' = {l+z)Eg{z) 



, , dEgiz) , 

y{z) = —l^ = {i + z)e^v 



Ho Jo + z' J^ + Qa V dE' J E'=ii+z')E,(z') 



(B5) 



(B6) 



where Eg{z) is an energy at epoch z. In case of the highest energies the intergrand in Eq. (jB6p can be expanded in 
powers of z', and after simple transformations one obtains 

dEg/dE^b{Eg)/b{E), 
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valid with 10% accuracy at _E > 4 x 10^^ eV. 
Eq. (HI]) coincides with Eq.^Se) from [U. 

We shah give now another derivation of Eq. (jB6P based on the exact sohition to the kinetic equation for propagation 
of UHE protons with continuous energy losses. This solution automatically includes dEg/dE term which coincides 
with that given by Eq. (|B6p . This equation for density of UHE protons, np{E, t), reads 



-npiE,t)-— b{E,t)np{E,t) ^Qg{E,t)=Q, 



(B7) 



with h{E, t) = EH{t) + b{E, t), where EH{t) describes adiabatic energy loss and b{E, t) - those due interaction with 
CMB. 

Written in the equivalent form 



^np{E,t) - biE,t)^np{E,t)-npiE,t)^^^^ - Qg{E,t) = 0, 



(B8) 



this equation can be solved with help of auxiliary characteristic equation dE/dt = —b{E,t), which solution is £{t) — 
Eg{Ei, ti,t), where like in Section [ll Bl Eg is the generation energy at age t, ii E — Ei at ti. When energy of a proton 
is taken on the characteristic E — £{t), the second term in the l.h.s. of Eq. (|B8p disappears and the solution found 
with help of integration factor is 



np{£,t) 



J' dt'Qg[£{t'),t']exp J'dt"-^b{£{t"),t") 



Introducing the red-shift z as a variable, dt/dz = —H ^(z)/(l + z), we obtain for present epoch (z = 0): 



np{£) 



dz' 



Using 

one finally gets 
np{£) = 



{l + z')H{z') 
db{£,z 



-Qg [£{z'),z']exp 



dz" db{£,z") 
(1 + z")H{z") d£ 



(B9) 



(BIO) 



d£ 



Hiz) + il + zY 



dboiE' 



dE' 



dz' 



{l + z')H{z') 



[£(z'),z'](l + ^')exp 



1 



dz"- 



E' = {l+z)£(z) 



ll\2 



{l + z") 



fdbo{E') 



^n^ii + z")^ + nA V dE' J E'^(i+z'.)£(z-') 



(Bll) 

Since £{z') — Eg{z'), one easily recognizes inside the integral the term dt'Qg{Eg,t') and dEg/dE as given by 
Eq. (|B6ll . 



APPENDIX C: KINETIC EQUATION SOLUTION VS CONTINUOUS LOSS APPROXIMATION: 

PHYSICAL INTERPRETATION 

We shall compare here the kinetic equation solution for proton spectra, nkin{E), with that obtained in continuous- 
energy-loss (CEL) approximation, ncont{E), and study the physical meaning of their difference. 

Additionally to the full kinetic equation we shall consider here also the stationary kinetic equation with the 
photopion production only. Such equation describes realistically the spectrum at > 1 x 10^° eV, where pion energy 
losses dominate, and effects of CMB evolution are negligible. However, we formally consider this equation at the lower 
energies, too, to study the difference between kinetic equation solution and CEL approximation. 

The general kinetic equation is reduced to the stationary one with photopion interaction only, putting in Eq. (1161) 
dn/dt = 0, H{t) = 0, T = To and bpair = 0. 

Then the stationary kinetic equation reads 

-PiE)np{E)+ J dE'P{E',E)np{E') + Qg,^{E) = Q, (CI) 

E 



35 



with P{E) and P{E' , E) given by Eqs. (|17p and ([18]). Introducing variable x = E/E', we rearrange the regeneration 
term into 

1 

^ j dx (^^^ P{E/x, x)np{E/x). (C2) 



Expanding the integrand into the Taylor series in powers of {1 — x) one observes the exact compensation of the zero 
power term with —P{E)n{E). The first power term results in the stationary continuous energy loss equation 

^ [b(E)np{E)] + Qge^iE) = 0, (C3) 



and account of (1 — x)^ term gives the stationary Fokker-Planck equation 

Qgen{E)=0, (C4) 



_d_ 



b{E)np{E) + A {E'D{E)nj,{E)) 



where the term with diffusion coefficient D{E) in momentum space describes fluctuations. 

This result reproduces the general proof (see e.g. [t^, p. 89) for the case when all processes go with a small energy 
transfer. 

Let us now consider the analytic solution to Eq. (|C1[) with -Emax — > oo in the asymptotic limit E 3> Eqzk- We shall 
use here {and only here) a toy model with simplified assumptions about interactions, namely with one-pion production 
p + 7 ^ TT + A^, with isotropic pion distribution in cm. system and with asymptotic cross-section ap~f — > const, when 
Ec —f oo. With these assumptions we readily obtain 

P{E) cGp^n^, P{E', E) cap^n^/E', (C5) 

and using Qgen{E) = QqE^^h ^ we arrive at 

nMn{E) = Q^E-^^. (C6) 

For continuous energy loss equation (|C3[) . using b{E') = (£" — E)c(jp^n^ and {E' — E) = \E' , valid for isotropic 
pion distribution in cm. system, we obtain 

ncont{E) ^ ^^QoE-<^ . (C7) 

Therefore, the asymptotic ratio is 

nkrnjE) ^ 7g 

n.ont{E) 2 ' ^"^""^ 

valid also for 7g = 2. This result has been obtained earlier in |l2Cll |. Of course, at the considered energies nothing but 
fluctuations are responsible for ratio (jCSp . This ratio, being asymptotic, can be used, however, as an estimate of role 
of fluctuations at i? > > Eqzk ■ 

Let us now come over to the case of finite Emax and consider E very close to E^ax ■ 

The regeneration term in Eq. (jCl|l tends to zero and 

nMn{E) ^ —— -, (C9) 

^ \-^raax } 

having non-zero value at -B = Emax ~ e and being zero at E = Emax + e, where e is infinitesimally small energy. 
For particle density aX E = Emax — e in continuous loss approximation one obtains from Eq. (|C3p 

1 



ncont{Emax) = 77T^ T / Q^E^^^ dE ^ > 0. (CIO) 



Therefore, nktn{E)/ncont{E) ^ oo as ^ ^„ 
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This result has a clear physical meaning. In our model sources are homogeneously distributed and from nearby 
sources the particles with energy E,nax can arrive at the observational point with the same energy due to fluctuations. 
But in CEL approximation they cannot do it, since they loose energy at any distance traversed, whatever small it is. 
The described effect influences the ratio r{E) = nkin{E) / ncont{E) at all energies close enough to Emax- 
Let us now come over to the solution of Eq. (jCip with interaction given by the detailed calculations as described in 
Section|Vl The calculated ratio r{E) = nkin{E) /ncont{E) is plotted in Fig. [23] for -Emax = 1 x 10^^ eV and for three 
values of 7g equal to 2.3, 2.7 and 2.9. One may observe in Fig. [23] some features described by our toy model: the 
ratio r{E) increases with 7^ as it should according to Eq. (|C8p and it tends to 00 at Emax- For spectra of interest for 
the present work with 2.3 < 7g < 2.8 the ratio r{E) at 1 x 10^^ < E <1 x 10^^ eV does not exceed 10 %. Note that 
in Fig. [15] it is smaller at low energies because the non-fluctuating pair-production process is included there. 

The ratio of the solutions for kinetic equation and for CEL equation includes fluctuations, but it includes also 
some other effects. The kinetic equation tends to CEL equation (|C3p in the limit of small 1 — x. But in fact for all 
reasonable energies 1 — a: is limited by minimum value 771^/(771,^ + 77ip) « 0.13, which cannot be considered as small 
value for the discussed accuracy. The same argument is valid for comparison of Monte-Carlo simulation with CEL 
approximation: why they should coincide if the energy transfer process is discrete? However, we do see from Fig. [53] 
that the ratio r{E) differs from 1 just by a few percent, when energy becomes less than 1 x 10^^ eV. 

Appendix [Dl gives an example when kin/CEL ratio at energy close to 1 x 10^" eV is explained not by fluctuations. 
At energy E k, Q x 10^^ eV the kin/CEL ratio is affected by an interplay between pair- and photopion production. 

We shall now discuss the comparison of our calculations with Monte-Carlo (MC) simulations. 

As was emphasized in Section [V] the kinetic equation solution must coincide exactly with MC simulation when the 
number of the MC tests tends to infinity, the distribution of sources is homogeneous, 7^ and £^max are the same and 
interactions are included in the identical way. D. De Marco and M. Kachelriess run their MC programs specially for 
comparison with our calculations using 7^ = 2.7, -Emax = 1 x 10^^ eV and homogeneous distribution of sources. They 
used SOPHIA interaction model which is not identical to our interaction model especially at the highest energies 
-E > 1 X 10^-'^ eV, but the average energy losses coincide well in both models (see Fig.[T]). The comparison of MC/CEL 
and kin/CEL ratios are given in the Table below. The data of MC are kindly provided by D. De Marco (D) and M. 
Kachelriess (K). 

From Table [III one can see the tendency of smaller kin/CEL ratios in comparison with MC/CEL, which can be due 
to different interaction models used in these calculations. We are planning the joint work with D. De Marco and M. 
Kachelriess on the detailed study of the discussed effects by methods of kinetic equation and Monte-Carlo. 

The main conclusion of our study is that fluctuations in photopion production modify only weakly the spectra with 
different 7^ at energies up to 1 x 10^^ eV (see Fig. 
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TABLE II: Comparison of kin/CEL and MC/CEL ratios 



E (eV) 


1 X 10^° 


1 X 10^^ 


1 X 10^^ 


kin/CEL 


1.05 


1.1 


1.3 


MC/CEL {D) 


1.16 


1.16 


1.25 


MC/CEL (K) 


1.1 


1.3 





APPENDIX D: THE SECOND DIP 



As explained qualitatively in Section FlII DI the second dip appears due to breaking of compensation between particle 
exit and regeneration given by the first and second terms in Eq. (jCip . This is a narrow feature near the energy -Ecq2 
which is produced because of the sharp (exponential) increase of photopion energy losses with energy, seen in Fig.[TJ At 
energy E just slightly below £^cq2 the continuous (adiabatic and pair-production) energy losses dominate, and spectrum 
is universal. At energy slightly above i?cq2 the continuous energy losses are small, and spectrum is determined by 
Eq. ()C1|) with the high accuracy compensation between the absorption, —P{E)np{E), and regeneration terms. At 
E Eccft the continuous energy losses break this compensation, increasing the absorption term, and triggering thus 
the appearance of the dip. We shall study here this triggering mechanism quantitatively. 

Coming back to the basic kinetic equation (|16p . we rearrange the terms connected with expansion of the universe 
(proportional to H{t)) and with pair-production energy losses, including them in Pgff{E,t). We obtain 



dnp{E,t) 
dt 



-Peff{E,t)np{E,t)+ J dE'P{E',E,t)np{E\t)+Qgen{E), (Dl) 



E 



where Peff{E,t) ^ P{E,t) + Pcont{E,t), with P{E,t) given as before by Eq. ^ and 

P,ont{E,t) ^ 2H{t) - [Pp,,r{E,t)+H{t)] ^^^^^ - (D2) 

As was described above, Pcont{E) breaks the compensation between —P{E)np{E) and the regeneration term in 
Eq. (|D1[) , triggering thus the modification of n^in {E) . It is convenient to introduce the auxiliary trigger function 
defined at t = to as 

TiE) - / PeffiE)/Pcont{E) at E<E, 
^ ' \ Peff{E)/P{E) E>E,, ^ ' 

with Ec = 6.1 X 10"'^^ eV, determined from equation P{E) = Pcont{E) and being equal to £'oq2. The trigger function 
describes how Peff{E) in Eq. (jPip is changing from Pcont{E) aX E <^ Ec, where T = 1, to P{E) at E Ec, where 
T = 1 as well. The calculated trigger function is plotted in the left panel of Fig. [Ml In the calculations we used for 
np{E) the universal spectrum, because as Fig. 1151 shows, its distortion is small. The trigger function T reaches 1 at 
i!^ < 3 X 10^^ eV, and therefore at these energies Pcff{E) w Pcont{E); the regeneration term is small at these energies, 
too, and thus nkin{E) w UcontiE) holds. At > 1 x 10^" eV T w 1, Pcff{E) w P{E) and nkin{E) becomes the 
solution of Eq. (|C1|) with photopion energy losses only and with nkin{E) /ucontiE) as shown in Fig. 1231 The maximum 
deflection of n{E) from these two regimes reached at Ec — 6.1 x 10^^ eV. For the steep spectrum n{E), which is the 
realistic case, the second term in Eq. (|D2p is positive and large in comparison with other terms, Pcont{E) > 0, the 
effective exit probability P{E) + Pcont{E) increases, and thus nkin{E) decreases. The triggering mechanism predicts 
that Ec does not depend on 7g and -Emax- 

The exact calculations for rikin{E) / ncont{E) are shown in the right panel of Fig. [Ml In fact, this is just the ratio 
of the modification factors from Fig. [T5I Cfluct.' and 'cont.loss.' curves), shown with a strong multiplication in order 
to demonstrate the fine structure of this function. One can observe the narrow dip at i52dip = 6.3 x 10^^ eV, which 
is only marginally seen in Fig. 1151 The second dip appears in the solution of the Fokker-Planck equation, too. 

The minimum of the dip coincides well with the maximum of the trigger function = 6.1 x 10^^ eV, and the small 
difference between them is caused by the fact that Ec is calculated for t — to, while for E2dip the contribution from 
epochs with t < t^ is present. The widths of both features are numerically similar. As the exact calculations show, 
the position of the dip does not depend on 7^ and ii^max, in accordance with triggering mechanism. 
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FIG. 24: In the left panel the trigger function T{E) is displayed as function of energy. In the right panel the ratio k{E) = 
nkin{E) /ncont{E) is plotted for three values of -Bmax- Proton density nkin(E) is found as a solution of the total kinetic 
equation (|16p . and ncont{E) is the universal spectrum calculated in the continuous energy loss approximation for homogeneous 
distribution of the sources. The position of the second dip, E — 6.3 x 10^^ eV, coincides well with the maximum of trigger 
function Ec = 6.1 x lO'^^ eV, and widths of both features are similar. 



We would like to emphasize the similarity of the first and the second dip. The first dip starts at energy i?eqi, 
where pair-production energy losses become equal to the adiabatic energy losses. The second dip occurs near energy 
Eeq2, where pair-production and pion-production energy losses are equal. Both dips are the faint features, not seen 
well when measured spectrum is exposed in a natural way JohsiE) vs E. The first dip is clearly seen in ratio of 
measured spectrum to unmodified spectrum, Johs{E) / J^^nmiE) , the second dip is predicted to be visible in the ratio 
of the measured spectrum to the smooth universal spectrum, Johs{E)/ Juniv{E). The latter is predicted as ratio 
{E)/ncont{E) shown in Fig.[Ml 



APPENDIX E: THE GENERATION FUNCTION FOR CASE OF DISTRIBUTION OF SOURCES OVER 

MAXIMAL ENERGIES 

We calculate here analytically the generation function Qgen{E), which gives the number of produced particles per 
unit comoving volume of universe and per unit time at z = 0, in case the sources are distributed over maximum of 
acceleration e = E^ax- 

We consider first the acceleration by non-relativistic shock, assuming the power-law spectrum of accelerated particles 
~ E~'^ with exponential suppression exp(— i?/e) at -E > e. Namely, we take the generation function for a single source, 
Qgen{E), which givcs the number of particles with energy E produced per unit time, in the form 

1sen{E,e) = -—^^ SV9en{E,e), (El) 

where Emin 1 GeV is the minimum energy in acceleration spectrum, Lp{e) is the luminosity of a source in the form 
of accelerated particles with energy E, and 



' e cxp(l — £^/e) at E > e. 



Then the generation function per unit volume Qg{E) is 



Q,en{E)=l ^^^de 2iT,f^ . V,en{E,e), (E3) 



where ris(e) = ^^[e, Lp{e)] is the space density of the sources. 
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We introduce also the spectral emissivity, = ns{e)Lp{e) and the total emissivity Cq = J C{e)de. 
Qgen{E) Can bc evaluated from Eqs. (|E3p and (|E2p as follows: 

for E < Ermn, 

Qgen{E) ~ E-^- (E4) 
for Emin < E < Emax the main contribution to Qgen{E) is given by 

r^^dECiE), (E5) 
Je 

while the integral from Emin to E results in 

QfliE)^E-'CiEy, (E6) 

and for E > Emax 

Qgen{E) - exp(-£;/e 

max )• (E7) 

The physically interesting regime is given by Emin < E < Emax- To reproduce Qgen{E) ^ E~'^-'^ needed for the best 
fit of the dip, £(e) from (jE5|) and (jE6P must be power-law function ^ E~^ with /3 = 1.7. 
In this case we obtain 

at E <^ Emin 

QgeniE) - { E'^'^ at £„„„ < £; < (E8) 

^t E > Emax-; 

which coincides with the phenomenological broken generation spectrum (|20p . 

The function QgeniE) from (|E8p can be easily normalized using the condition J EQgen{E)dE = Cq. 
One may generalize this derivation for the case 

^ \ at E<E 

^-"^"^'^)^|^;-(-")/(f) at ^;>e, 

where f{x) is a function which diminishes with x like or faster and /(I) — 1. 
Similar to ([El]) and (jE3| we have 

QgeniE, e) = aE^-^i„Lp{E)ipgeniE,E), (ElO) 
Qge«(£',e) = aE"^i„ J dEC{E)ipgen{E,E). (Ell) 

Using the properties of /(x) we can prove again that the power-law function Qgen{E) ~ E^'^n in the interval 
Smin < E < Emax uccds the powcr-law function £(e) ^ e^^. 
Then using again the properties of /(x) we obtain 




-vain 



£:-(2+a) E <En 

i{E) ^ { at £„„;„ <E< Emax (E12) 

I i.E I Emax) '^jt -£/ > ^maX7 

We are interested in particular in case of acceleration at relativistic shock, when a « 0.2 — 0.3. The case of 
7g = 1 + a + /3 = 2.7 results then in /3 = 1.7 - a w 1.4 - 1.5. 
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